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NUMERICAL  CALCULATION  OF  SURFACE  WAVE  EFFECTS 


ON  MARINE  STRUCTURES 


INTRODUCTION 

The  effects  of  ocean  waves  on  man-made  structures  are  of  great  practical  interest,  but  until  quite 
recently  the  complexities  of  the  problems  have  limited  theoretical  calculations  to  idealized  flows  and 
structural  geometries.  In  the  absence  of  flow  separation  the  complexity  arises  primarily  from  the  pres¬ 
ence  of  the  free  surface,  particularly  if  wave-breaking,  surging  or  slamming  are  involved.  Not  only  does 
a  regular  free  surface  alone  easily  dominate  the  flow,  but  also  it  usually  is  accompanied  by  directional 
and  irregular  waves,  and  wave/current  interactions.  The  development  of  means  to  calculate  the  effects 
of  these  complex  flows  poses  intractible  problems  not  only  analytically,  but  also  numerically,  and  only 
the  latest  generation  of  computers  has  made  practical  calculations  even  remotely  feasible.  This  is  par¬ 
ticularly  true  for  parametric  studies  which  require  the  calculation  of  time-dependent  flow  fields  and  sur¬ 
face  pressure  distributions  for  a  variety  of  body  shapes  over  a  spectrum  of  wave  heights,  wavelengths, 
and  depth  conditions. 

A  research  program  is  underway  at  the  Naval  Research  Laboratory  to  develop  means  for  deter¬ 
mining  the  effects  of  surface  waves  on  submerged  and  partially  submerged  structures  in  the  ocean.  One 
principal  thrust  of  the  program  has  been  to  develop  a  capability  for  the  numerical  calculation  of  non¬ 
linear  free-surface  waves,  with  verification  through  comparisons  with  existing  wave  theories  and  force 
models,  and  with  experiments  in  a  wave  channel.  The  initial  stages  of  this  aspect  of  the  research  pro¬ 
gram  have  centered  upon  extending  the  capabilities  of  the  computer  code  SPLISH  which  was  developed 
at  NRL  by  Fritts  and  Boris  (1975,  1977b). 
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While  there  are  many  approaches  that  might  be  taken  for  the  numerical  calculation  of  free-surface 
waves,  e.g.  the  finite-element  work  of  Bai  and  Yeung  (1974)  and  Wellford  (1978),  and  the  marker- 
and-cell  codes  of  Nichols  and  Hirt  (1977),  there  are  distinct  advantages  in  using  a  technique  which 
naturally  tracks  the  wave  flow  field  and  which  employs  a  grid  by  which  the  free-surface  and  structural 
boundaries  can  be  readily  defined.  The  non-orthogonal  curvilinear  coordinate  systems  developed  by 
Shanks  and  Thompson  (1977)  can  accommodate  irregular  boundaries  without  interpolation  but  require 
transformation  calculations  at  each  time  step.  SPLISH  couples  the  advantages  of  a  Lagrangian  approach 
to  a  triangular  grid  which  can  be  restructured  easily  in  a  conservative,  reversible  manner  to  meet  the 
needs  of  calculations  of  transient  phenomena. 

This  report  presents  a  description  of  the  numerical  method  upon  which  SPLISH  is  based.  The 
numerical  results  are  compared  with  linear  and  fifth-order  wave  theory  results  and  the  numerical  and 
theoretical  results  are  compared  with  experimental  data.  For  one  range  of  conditions  the  agreement  is 
good  between  the  present  numerical  results  and  results  from  both  theory  and  experiment,  and  this 
agreement  seems  to  validate  the  SPLISH  computer  code.  However,  the  results  obtained  for  another  set 
of  conditions  characterized  by  large  wave  reflection  illustrate  the  need  for  alternate  boundary  conditions 
in  the  code. 

THE  NUMERICAL  METHOD 

Ideally,  calculation  of  the  effects  of  ocean  waves  on  a  structure  should  incorporate  boundary  con¬ 
ditions  which  permit  a  spectrum  of  wavelengths  and  an  adequate  range  of  wave  heights  for  different  sea 
conditions.  The  modeling  of  transient  phenomena  such  as  surges,  wave-breaking,  and  sp>  ay  production 
is  necessary  in  the  long  term  since  these  phenomena  account  for  a  major  portion  of  wave  energy  dissi¬ 
pation.  The  SPLISH  code  is  structured  to  permit  calculations  in  the  presence  of  these  complex  effects 


under  certain  conditions. 
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For  two-dimensional,  incompressible,  and  inviscid  fluid  flow  —  which  is  the  situation  of  concern 

for  calculating  free-surface  wave  flows  in  the  ocean  —  the  governing  equations  are  almost  deceptively 

simple.  The  equations  are  the  momentum  equation 

d\  „ 

p—  ~-Vp  -  pgy  (1) 

and  the  continuity  equation 

V  V  -  0.  (2) 

Closure  of  the  system  of  equations  is  given  by  the  boundary  conditions.  At  solid  boundaries,  such  as 
the  bottom  of  the  channel,  the  normal  velocity  is  given  by 

Vk  -  0.  (3a) 

At  the  free  surface  the  pressure  (relative  to  atmospheric)  is 

P  =  0.  (3b) 

The  periodic  boundary  conditions  used  in  the  current  version  of  SPLISH  require  (at  the  sides  of  the 

computational  region) 

/(0)  -  f(k).  (4) 

where  /  is  an  arbitrary  dependent  variable  and  X  is  the  wave  length,  which  limits  consideration  to 

domains  of  integral  numbers  of  wavelengths.  At  the  free  surface  y  —  tj,  dynamic  and  kinematic  boun¬ 
dary  conditions  on  the  wave  flow  are  derived  by  satisfying  Bernoulli’s  equation  and  by  requiring  that  the 
substantial  derivative 

~  (y  -  v)  =  0.  (5) 

Here  ij  is  the  time-dependent  elevation  of  a  point  on  the  free  surface  (see  Appendix  B). 

The  governing  equations  are  solved  by  a  finite-difference  method  which  was  developed  especially 
for  use  on  a  triangular  grid  which  moves  with  the  flow.  The  finite-difference  method,  the  Lagrangian 
approach,  and  the  triangular  grid  are  discussed  in  depth  in  Appendix  A  and  are  summarized  below. 

A  Lagrangian  grid  will  distort  in  any  non-trivial  flow  field,  and  as  grid  distortion  becomes  severe 
the  calculation  quickly  loses  accuracy.  However  a  triangular  grid  can  be  manipulated  locally  in  several 
ways  to  extend  realistic  calculations  of  transient  flows.  Each  grid  line  represents  a  quadrilateral  diago¬ 
nal,  and  the  opposite  diagonal  can  be  chosen  whenever  vertices  move  in  the  flow  to  positions  which 
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favor  that  connection.  Such  a  reconnection  involves  just  the  four  vertices  describing  the  quadrilateral. 
No  fluid  moves  relative  to  the  quadrilateral,  eliminating  one  form  of  numerical  dissipation.  Vertices 
may  also  be  added  or  deleted  to  preserve  the  desired  resolution  by  local  algorithms  which  involve  only 
those  vertices  in  the  vicinity  of  the  grid  anomaly.  Major  advantages  of  this  technique  are  that  the  algo¬ 
rithms  can  be  conservative,  they  permit  a  minimum  of  numerical  dissipation  and  yet  they  require  very 
little  computer  time  since  most  of  the  grid  remains  unaltered. 

A  basic  technique  of  differencing  the  equations  over  such  a  grid  for  the  case  of  inviscid 
incompressible  flow  is  given  in  Appendix  A.  The  current  version  of  the  code  employs  centered  implicit 
integration  to  advance  the  physical  variables  in  time  and  has  been  shown  to  be  second-order  accurate. 
Voriicity  is  conserved  identically  and  the  pressure  is  advanced  to  conserve  the  divergence  about  each 
vertex.  These  new  pressures  are  then  used  to  advance  the  velocities  and  update  the  grid  positions. 

Although  much  of  the  work  thus  far  has  centered  on  the  finite-difference  techniques  and  on  the 
grid  restructuring  algorithms,  the  problems  of  realistic  boundary  conditions  and  initial  conditions 
remain.  Inflow  and  outflow  boundary  conditions  which  do  not  adversely  affect  the  rest  of  the  computa¬ 
tional  region  over  a  spectrum  of  wavelengths  are  goals  which  have  not  yet  been  reached  even  on  more 
standard  grids.  Similarly,  initialization  procedures  for  nonlinear,  transient  flows  involving  a  free-surface 
have  been  investigated  only  minimally  and  on  a  case-by-case  basis. 

Since  techniques  for  boundary  conditions  and  initialization  procedures  of  free-surface  waves  are 
still  in  an  early  stage  of  development,  the  present  results  necessarily  are  restricted  in  their  range  of 
application.  The  choice  of  periodic  boundary  conditions  restricts  the  calculation  domain  to  an  integral 
number  of  wavelengths  and  requires  that  the  domain  is  one  element  of  an  array  of  identical  domains. 
For  many  conditions  of  wave  flow  this  yields  a  quite  realistic  situation  and  there  is  no  concern  for 
numerical  damping  at  radiative  boundaries.  This  restriction  on  the  problem  has  in  fact  permitted  con¬ 
siderable  progress  to  be  made  in  the  numerical  calculation  of  free-surface,  progressive  waves,  even  for 
cases  with  obstacles  such  as  bottom-seated  half-cylinders  in  the  flow.  The  present  findings  show  that 
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when  wave  reflection  from  an  obstacle  is  small,  the  periodicity  of  domains  does  not  prevent  good  agree¬ 
ment  between  numerical  results,  experimental  data  and  theoretical  results.  In  fact,  the  agreement  is 

quite  good.  This  is  an  unexpected  result,  since  the  blockage  due  to  the  obstacle  can  be  quite  large  and 
deviations  from  purely  periodic  flow  would  be  expected. 

However,  if  wave  reflection  from  an  obstacle  is  significant,  the  domains  to  the  left  and  right  of  the 
center  computational  domain  have  a  very  large  influence  on  the  numerical  results.  For  the  high 
reflection  case  considered  here,  the  periodicity  of  domains  leads  to  a  solution  of  a  problem  which  is 
different  from  the  problem  of  interest.  This  is  not  surprising  since  the  problem  posed  numerically  is 
quite  different.  Even  so,  some  of  the  numerical  results  appear  to  be  useable.  The  high  reflection  case 
does  emphasize  the  desirability  of  future  research  to  include  the  formulation  of  continuative  (i.e.,  non¬ 
periodic)  boundary  conditions  which  would  permit  more  realistic  calculation  of  the  wave  flow  over  a 
single  obstacle. 

The  calculation  is  started  by  a  sinusoidal  pressure  pulse  on  the  free  surface  which  initiates  a  stand¬ 
ing  wave.  At  the  quarter  period  of  the  standing  wave  ( t  —  7/4  sec)  a  second  pulse  (phase  shifted  a 
quarter  wave  length)  is  applied  to  form  the  progressive  wave.  This  initialization  procedure  may  not  be 
the  most  appropriate,  especially  when  an  obstacle  is  in  the  flow  field,  but  it  is  simple,  easy  to  apply,  and 
exhibits  good  numerical  stability  for  standing  waves  (Fritts  and  Boris,  1977a)  and  for  progressive  waves 
in  a  straight  channel.  Also,  it  allows  the  periodic  boundary  conditions  to  be  retained  and  introduces 
minimal  grid  distortion  and  variation  in  the  triangle  areas. 

RESULTS  AND  DISCUSSION 

In  the  following  sections  we  have  three  primary  objectives.  The  first  objective  is  to  validate  the 
numerical  method,  i.e.,  the  computer  code  SPLISH,  for  calculating  progressive  free-surface  waves.  The 
second  objective  is  to  demonstrate  the -applicability  of  the  numerical  results.  The  third  objective  is  to 
demonstrate  the  potential  of  the  numerical  approach,  especially  for  those  cases  where  theory  is  limited 
by  approximations. 
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The  results  are  presented  in  the  following  two  sections.  First,  the  numerical  results  are  discussed 
and  compared  with  results  from  classical  wave  theory.  In  the  subsequent  section  the  present  numerical 
results  from  SPL1SH  are  compared  with  experimental  data  obtained  in  the  NRL  wave  channel.  Good 
agreement  is  obtained  between  the  numerical  results  and  various  theories  and  between  the  numerical 
results  and  the  experiments  when  the  numerical  simulation  appropriately  models  the  theory  or  experi¬ 
ment. 

Numerical  Results  and  Comparison  With  Theory 

Results  from  calculations  for  three  sets  of  conditions  are  presented  and  a  comparison  is  made  with 
results  obtained  from  wave  theories.  Although  the  numerical  formulation  is  not  limited  to  linear  or 
nearly-linear  waves,  the  calculations  were  performed  for  small  amplitude  waves  to  permit  meaningful 
comparison  with  theory. 

The  first  conditions  considered  are  for  a  two-dimensional  wave  channel  with  no  obstruction  to  the 
flow.  For  single  frequency  waves,  the  periodic  boundary  conditions  are  quite  appropriate  and  the  good 
agreement  with  the  theory  establishes,  in  part,  the  validity  of  the  numerical  calculations.  The  next  con¬ 
ditions  considered  are  for  a  half-cylinder  mounted  on  the  bottom  of  the  channel  with  wavelength  (X), 
depth  (d)  and  cylinder  radius  ( a )  chosen  to  yield  minimal  wave  reflection  from  the  half  cylinder, 
according  to  the  theory  of  Naftzger  and  Chakrabarti  (1979).  In  the  third  case,  the  values  of  X,  d and  a 
are  changed  so  that  the  amount  of  wave  reflection  is  significant.  The  numerical  results  for  this  case 
clearly  show  the  need  for  continuative  boundary  conditions  since  the  numerical  simulation  does  not 
properly  model  the  same  physical  situation. 

Uniform  Depth  Channel 

For  this  simplest  of  cases  to  consider,  the  computational  domain  is  divided  into  triangles  as  shown 
in  Figure  1.  Twenty-six  vertices  were  used  in  the  horizontal  direction  and  six  in  the  vertical  direction, 
for  a  total  of  156  vertices.  The  lack  of  direct  correspondence  between  vertices  and  triangles  is  discussed 
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in  Appendix  A  and  is  illustrated  here  as  2S0  triangles  form  the  computational  domain.  Once  the 
domain  is  divided  into  triangles,  the  calculation  is  started  by  a  sinusoidal  pressure  pulse  on  the  free  sur¬ 
face  which  initiates  a  standing  wave.  Then  at  the  quarter  period  of  the  standing  wave  (/  -  7/4  sec)  a 
second  pulse  (phase  shifted  a  quarter  wave  length)  is  applied  to  initiate  the  progressive  wave.  The 
advantages  of  this  initialization  procedure  are  simplicity,  lack  of  grid  distortion,  and  good  results  which 
have  been  obtained.  A  disadvantage  is  the  transients  from  the  pressure  pulses,  but  these  transients  do 
not  seem  to  cause  any  serious  problems.  The  transients  can  be  reduced  by  an  alternative  initialization 
procedure  wherein  the  free  surface  is  displaced  to  a  modified  sine  wave  (Penny  and  Price,  1952)  and 
the  progressive  wave  then  formed  by  a  single  pressure  pulse  at  /  =  0.  With  this  procedure  the  pressure 
transients  are  reduced  but  some  grid  distortion  is  introduced. 

The  triangular  grid  (for  depth  d  “  1  m  and  wavelength  X= 2.5  m)  plotted  for  three  points  in  time 
in  Figure  2  illustrates  a  potential  disadvantage  of  the  Lagrangian  approach.  The  grid  is  shown  at  times  t 
—  1.2,  1.6  and  2.0  sec  (or  at  time  steps  31,  41  and  51)  and  the  wave  motion  from  left  to  right  is  clearly 
shown.  Since  this  particular  calculation  was  performed  without  grid  reconnections,  the  movement  of 
fluid  particles  (i.e.  triangle  vertices)  with  the  progressive  wave  introduces  some  very  obvious  distortion 
into  the  grid.  With  further  time  increments,  the  grid  distortion  would  become  even  more  pronounced 
and  surface  vertices  could  be  carried  far  to  the  right  of  vertices  initially  below  them.  Since  pressure 
gradients  are  calculated  over  triangles,  it  is  necessary  that  a  subsurface  vertex  be  connected  to  surface 
vertices  near  it  and  not  to  vertices  which  the  flow  has  carried  far  to  the  right.  As  discussed  in  Appen¬ 
dix  A,  the  reconnection  algorithms  reestablish  the  grid  connections  so  that  nearby  vertices  are  con¬ 
nected.  Triangles  which  are  carried  off  to  the  right  by  the  flow  are  replaced  at  the  left  by  identical 
triangles  under  the  periodic  boundary  conditions.  Such  grid  reconnections  and  triangle  movements 
eliminate  most  of  the  grid  distortion.  Figure  3  shows  the  grid  at  t  ~  2.0  sec  for  a  similar  calculation 
with  grid  reconnection  permitted  and  with  most  grid  distortion  eliminated.  The  reconnection  algo¬ 
rithms  are  crucial  to  the  success  of  a  Lagrangian  procedure,  especially  when  there  is  significant  long 
term,  net  motion  of  some  fluid  particles  relative  to  other  particles.  The  simplification  achieved  with  the 
periodic  boundary  conditions  is  also  illustrated. 
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The  calculations  with  SPL1SH  for  the  uniform  depth  case  also  gave  good  agreement  with  classical 
wave  theory  for  the  wave  period  and  for  the  mean  drift.  Previous  calculations  with  SPLISH  (Fritts  and 
Boris,  1977)  had  given  wave  period  values  for  a  standing  wave  which  agreed  very  well  with  theory.  In 
fact  the  period  converged  to  the  theoretical  value  as  the  grid  step  size  was  reduced.  Two  cases  for  a 
progressive  wave  (depth  d  =  1  m  and  wavelength  X  =2.5  m  and  thus  d/k  —  0.4)  were  considered  and 
good  agreement  was  obtained  with  the  theory  for  both  the  wave  period  and  the  mean  drift  motion.  For 
these  values  of  d  and  X  the  dispersion  relation  in  Appendix  B  (Equation  B6)  yields  a  wave  period  T  — 
1.28  sec.  For  deep  water  waves  the  period  would  be  T  =  1.27  sec,  so  for  this  case  the  progressive 
waves  are  very  nearly  deep  water  waves  and  the  deep  water  approximation  can  be  used. 

The  path  over  which  a  surface  particle  moves  during  several  wave  cycles  is  plotted  in  Figure  4. 
The  wave  amplitude  is  H  =  0.047  m,  or  H/k  —  0.019,  so  that  the  linear  wave  potential  given  by  Equa¬ 
tion  B5  is  applicable.  The  SPLISH-generated  path  of  a  surface  particle  is  shown  in  the  figure  by  the 
solid  line,  and  the  particle  path  predicted  by  integrating  Equations  B12  and  B13  is  shown  by  the  dotted 
line.  (Note  that  the  horizontal  scale  is  stretched  relative  to  the  vertical  scale.)  The  results  obtained  by 
the  two  methods  are  generally  in  good  agreement,  but  the  linear  theory  predicts  slightly  greater  mean 
drift  over  three  cycles  of  the  wave.  The  predicted  drift  of  the  surface  particle  is 
A  3c,  =  0.09  m  (At  =  3D  from  the  computer  code,  whereas  Ax;  =  0.105  m  (A t  —  3D  from  Equation 
B17.  The  average  calculated  wave  period  is  Dplish  =  132  sec  over  three  cycles  and  this  is  in  good 
agreement  with  the  wave  period  obtained  from  the  dispersion  relation  7Vt0kes  =  1 .28  sec.  As  in  the 
standing  wave  case,  Dplish  should  converge  to  Tstokes  for  decreasing  grid  size. 

The  Lagrangian  particle  paths  for  vertices  initially  at  y  *=  0  and  y  —  -0.2  m  are  plotted  in  Figure  5 
for  a  wave  amplitude  H  =  0.066  m  (  H/k  -  0.026  )  at  the  surface.  For  this  series  of  computations 
the  wave  period  is  rSPL|SH  “  134  sec  averaged  over  the  five  cycles  in  the  figure.  The  corresponding 
wave  period  from  Equation  B6  is  T  —  1.29  sec,  again  for  d/k  —  0.4.  The  vertical  particle  displacement 
H  at  a  depth  of  y  —  -0.2  m  (-0.08X)  is  H'  —  0.042  m,  which  results  in  the  ratio  H'/H  ™  0.64.  From 
the  linear  theory  the  particle  displacement  H'  at  the  depth  y  is  given  by 
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=  J^sinh  k(y  +  d) 
sinh  kd 

Applying  this  equation  for  ky  =  —0.5  and  kd  =  2.51  as  in  Figure  5  results  in  the  ratio  H'/H  **=  0.60. 
which  is  in  good  agreement  with  the  SPLISH-generated  particle  displacements  over  five  wave  cycles  at  y 
=  —0.2  m. 

The  mean  drift  of  the  fluid  particles  in  the  direction  of  the  wave  motion,  computed  from  Equation 
B17,  is  plotted  in  Figure  5  and  on  the  average  agrees  well  both  on  and  below  the  surface  with  the 
SPLISH-generated  particle  paths  shown  there.  While  the  results  at  .y  =  -0.2  m  agree  on  the  average 
over  five  wave  periods,  it  should  be  noted  that  the  particle  drift  per  cycle  varies  from  Ax,  =  0.014  m  to 
0.038  m  about  a  mean  value  of  A3c,  =  0.026  m  in  the  case  of  the  computer-generated  Lagrangian 
particle  paths  for  the  individual  cycles. 

These  results  establish,  at  least  in  part,  the  validity  of  the  numerical  approach  even  though  these 
particular  calculations  were  made  with  one  of  the  earlier  versions  of  the  SPLISH  code.  Subsequent 
developments  in  the  code  have  improved  the  numerical  results,  especially  the  uniformity  of  the  particle 
paths,  as  will  be  seen  below.  While  the  numerical  simulation  of  the  uniform  depth  channel 
configuration  closely  corresponds  to  a  realistic  physical  situation,  other  configurations  are  of  more  prac¬ 
tical  interest.  Two  such  configurations  are  discussed  next. 

Low  Reflection  Wave  Flow  over  a  Half-Cylinder 

One  configuration  which  is  of  some  practical  interest  is  the  wave  flow  over  a  bottom  seated  half¬ 
cylinder,  which  might  correspond  to  a  half  buried  pipe  or  a  storage  tank  on  the  ocean  floo'r.  This  case  is 

.* 

*-  ? 

also  relatively  easy  to  model  numerically.  However,  with  periodic  boundary  conditions  the  numerical 
simulation  corresponds  to  an  array  of  pipes  which  are  spaced  one  or  more  wavelengths  apart.  For¬ 
tunately,  there  are  combinations  of  wavelength,  depth,  and  cylinder  radius  for  which  only  a  very  small 
amount  of  the  incident  wave  is  reflected  by  the  half-cylinder.  For  example,  for  X  —  2.5  m,  d  —  1  m. 
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and  the  cylinder  radius  a  *  0,5  m  together  with  the  assumption  of  linear  theory,  Naftzger  and  Chakra- 
barti  (1979)  find  that  the  reflection  coefficient  R  —  0.03.  With  a  reflected  wave  in  this  case  having  only 
3%  of  the  amplitude  of  the  incident  wave  (and  thus  less  than  1%  of  the  energy),  neighboring  half¬ 
cylinders  should  have  negligible  effect  on  the  flow. 

A  series  of  calculations  was  made  for  the  low  reflection  conditions  described  above.  The  compu¬ 
tational  grid  for  one  such  case  (without  grid  reconnections)  is  shown  in  Figure  6.  Here  the  presence  of 
the  half-cylinder  is  indicated  by  the  upward  displacement  of  the  triangle  vertices  above  it.  In  subse¬ 
quent  calculations  the  amount  of  grid  distortion  was  reduced  by  permitting  grid  reconnections  and  by 
using  more  layers  of  vertices.  For  most  wave  amplitudes,  the  best  results  for  this  case  were  obtained 
when  9  rows  of  26  vertices  were  used  to  form  the  triangular  grid. 

Figure  7  shows  the  free  surface  contour  and  the  resultant  bottom  pressure  Pb  at  one  instant  in  the 
passage  of  a  wave  (from  left  to  right)  over  a  half-cylinder  (radius  a  =  0.5  m)  in  a  channel  (depth  d  — 

1  m)  for  a  wave  length  A  =  2.5  m.  The  finite  amplitude  of  the  wave  (//  *  0.045  m,  from  the  undis¬ 
turbed  free  surface)  causes  only  a  very  small  deviation  from  a  sinusoidal  surface  shape  in  spite  of  the 
presence  of  the  bottom  seated  half-cylinder.  Note  that  the  pressure  scale  increases  downward.  Since 
the  calculation  was  started  from  still  water,  the  pressures  at  step  1  correspond  to  the  hydrostatic  pres¬ 
sure  on  the  channel  bottom  and  on  the  surface  of  the  half-cylinder.  At  step  25  the  wave  gives  an 
increased  pressure  on  the  left  side  of  the  half-cylinder  and  a  reduced  pressure  on  the  right  side.  This, 

of  course,  is  responsible  for  a  net  force  (at  this  instant)  from  left  to  right.  At  other  steps  in  the 
calculation,  the  pressure  fluctuation  accurately  follows  the  passage  of  peaks  and  troughs  of  the  wave  as 

will  be  shown. 

The  full  wave  shape  and  the  bottom  pressure  distributions  are  well  documented  in  Figure  7,  but 
the  details  of  the  pressure  around  the  half-cylinder  are  not  shown  as  well  as  might  be  the  case.  In  Fig¬ 
ures  8-12,  the  domain  from  x  -  0.3  m  to  x  —  2.2  m  is  shown  for  times  of  t  —  0,  0.96,  1.28,  1.64  and 
1.92  sec.  As  in  Figure  7  the  scale  for  the  pressure  increases  downward.  The  scale  for  the  bottom  pres¬ 
sure  has  been  selected  so  that  at  t  *  0  the  Pb  curve  coincides  with  the  plot  of  the  bottom  contour.  The 
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calculation  was  started  at  still  water  level  (as  seen  in  Figure  8)  and  the  wave  travels  from  left  to  right. 
Figures  9-12  are  approximately  a  quarter  period  apart  and  the  progress  of  the  wave  can  be  seen  clearly. 

Figures  9-12  also  show  a  depth  dependence  of  the  pressure  fluctuation.  As  the  crest  and  the 
trough  of  the  wave  pass,  the  pressure  fluctuation  on  top  of  the  half  cylinder  (at  a  depth  of  O.S  m)  is 
approximately  twice  that  on  the  channel  bottom  (at  a  depth  of  1.0  m).  At  intermediate  depths  (i.e. 
along  the  sides  of  the  cylinder)  the  magnitude  of  the  pressure  fluctuation  is  between  that  on  the  chan¬ 
nel  bottom  and  that  on  the  top  of  the  cylinder.  These  figures  also  show  how  the  pressure  field  on  the 
channel  bottom  and  on  the  half  cylinder  correspond  to  the  passage  of  the  wave  on  the  free  surface.  In 
Figures  9-12  the  scale  for  the  pressure  curve  was  chosen,  as  mentioned  above,  so  that  at  t  =  0  the 
curve  for  Pb  coincided  with  the  contour  of  the  channel  bottom.  This  scale  suppresses  the  apparent 
magnitude  of  the  pressure  fluctuation.  The  fluctuation  in  Pb  can  be  accentuated  by  normalization. 

In  Figures  13  and  14  the  fluctuation  in  Pb ,  i.e.  AP,  is  normalized  by  its  maximum  value  (AFref) 
which  occurs  at  the  top  of  the  cylinder  as  the  wave  crest  passes  by.  Figure  13  shows  the  pressure  distri¬ 
bution  around  the  half-cylinder  at  times  t  *■>  1.28  and  t  **  1.92  sec  which  are  the  instants  when  the 
wave  crest  and  trough  pass  over  the  half-cylinder.  Figure  14  shows  the  pressure  distributions  around 
the  half-cylinder  at  t  =  0.96  sec  and  at  t  =  1.64  sec  which  correspond  to  times  when  the  wave  crest  is 
first  to  the  left  of  the  cylinder  and  then  to  the  right  of  the  cylinder.  In  these  two  figures  the  values  of 
A  P  from  the  SPLISH  numerical  calculations  are  shown  by  the  open  symbols.  The  solid  and  dashed 
curves  shown  in  these  figure  are  the  values  of  A  P  obtained  for  these  conditions  using  the  approach  of 
Chakrabarti  and  Naftzger  (1974)  which  was  based  on  Stokes’  fifth-order  wave  theory  and  the  assump¬ 
tion  that  the  effect  of  the  free-surface  on  the  reflected  wave  potential  could  be  neglected.  For  the 
present  low-reflection  conditions  (with  a  reflection  coefficient  of  3%)  such  an  assumption  is  quite 
appropriate.  For  these  conditions  the  SPLISH  results  are  in  good  agreement  with  the  fifth-order  theory 
results. 


We  might  note  here  that  for  this  case  there  is  properly  some  departure  from  the  results  that 
would  be  obtained  with  linear  wave  theory.  In  particular,  in  linear  wave  theory  the  magnitude  of  the 
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pressure  fluctuation  on  top  of  the  cylinder  (i.e.  at  0  =*=  0°)  is  exactly  the  same  when  the  wave  Crest  and 
the  trough  pass  by.  Such  is  not  the  case  with  higher  order  theories.  When  the  wave  is  described  by 
higher  order  theories,  the  crest  is  more  sharply  peaked  and  the  trough  is  more  shallow  than  for  the 
linear  theory.  Consequently,  the  distance  from  the  still  water  level  to  the  crest  is  greater  than  the  dis¬ 
tance  from  the  still  water  level  to  the  trough,  whereas  in  the  linear  theory  these  two  distances  are  equal. 
Thus,  for  A  P  =  0  corresponding  to  the  still  water  level,  the  magnitude  of  A  P  on  the  cylinder  for  the 
crest’s  passage  is  greater  than  for  the  trough’s  passage.  As  may  be  seen  in  Figure  13,  AP/APREF  is  1.0 
for  the  passage  of  the  crest  but  is  only  0.9  for  the  passage  of  the  trough.  The  results  of  SPLISH  and  the 
fifth-order  theory  agree  well. 

The  wave-induced  force  or  force  components  on  the  half-cylinder  can  be  obtained  by  integrating 
the  pressure  distribution,  such  as  those  in  Figures  13  and  14,  around  the  half-cylinder.  This  integration 
has  been  done  in  this  case  for  the  SPLISH  calculations  and  the  results  of  this  operation  are  shown  in 
Figures  15  and  16.  The  curves  for  the  force  components  are  started  at  t  —  0.76  sec,  after  the  initial 
transients  caused  by  the  pressure  pulses  which  form  the  progressive  wave  have  died  out.  The  first  two 
extrema  of  the  Fx  curve  in  Figure  15  correspond  closely  to  the  times  of  Figures  9  and  11  when  the 
crest  was  to  the  left  and  then  to  the  right  of  the  half  cylinder  and  with  the  free  surface  at  the  still  water 
line  above  the  half-cylinder.  Also,  the  first  two  extrema  of  the  Fy  curve  in  Figure  16  correspond 
closely  to  the  times  of  Figure  10  and  12  when  the  wave  crest  and  then  the  wave  trough  passed  above 
the  half-cylinder.  Also  shown  in  Figures  15  and  16  are  the  linear  theory  limits  for  the  maximum  values 
of  Fx  and  Fy  which  agree  well  with  the  SPLISH  results.  The  limit  values  were  calculated  using  a  wave 
amplitude  H  =  0.038  m  in  the  linear  theory  results  obtained  by  Naftzger  and  Chakrabarti  (1979).  The 
value  of  H  —  0.038  m  was  chosen  here  to  match  the  average  amplitude  (over  the  first  four  wave 
cycles)  of  the  wave  in  the  numerical  calculation.  This  was  necessary  since  the  wave  amplitude  (and 
thus  the  pressure  and  the  force  component  amplitudes)  show  a  decrease  with  time.  This  decrease  in 
wave  amplitude  (and  thus  wave  energy)  is  not  due  to  viscous  effects  in  the  governing  equations,  but  it 
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is  caused  by  an  effective  dissipation  of  energy  due  to  the  incomplete  convergence  of  the  successive- 
over-relaxation  (SOR)  method  used  for  the  solution  of  the  Poisson  equation  for  the  pressure.  The 
effect  of  this  decrease  in  wave  amplitude  will  also  be  seen  below. 

The  time  history  for  the  pressure  (in  kiloPascals  or  kN/m2)  at  a  vertex  on  the  top  of  the  half¬ 
cylinder  is  shown  in  Figure  17.  The  solid  curve  shows  the  SPLISH  data  and  the  open  symbols  are 
values  obtained  using  full  fifth-order  wave  theory  in  results  obtained  by  Chakrabarti  and  Naftzger 
(1974)  for  wave  flow  over  a  bottom  seated  half-cy Under.  The  pressure  here  represents  the  actual  gauge 
pressure  that  might  be  measured.  The  wave  amplitude  was  assumed  to  be  H  —  0.038  m  for  the  fifth- 
order  calculations;  this  had  been  done  in  calculating  the  linear  theory  limits  for  the  force  components 
plotted  in  Figures  15  and  16.  Here  the  agreement  between  the  numerical  calculations  and  the  fifth- 
order  calculation  is  quite  good.  With  the  finite  grid  size  employed  here  the  period  given  by  the  numeri¬ 
cal  results  is  about  4%  greater  than  that  from  the  theoretical  dispersion  relation,  so  that  the  theoretical 
and  numerical  results  shift  slightly  relative  to  each  other  with  time.  From  fifth-order  theory  T  =  1.269 
sec  and  from  linear  theory  T  —  1.274,  while  SPLISH  gives  here  a  value  of  T  =  1.30  sec  based  on  the 
intervals  between  maxima  or  between  minima  in  the  pressure-time  history  in  Figure  17. 

The  forces  and  pressures  on  the  half-cylinder  are  the  principal  results  of  the  numerical  calcula¬ 
tions,  but  at  times  the  necessary  information  about  the  quality  of  the  numerical  solution  is  not  given  by 
the  pressure  and  force  data.  Plots  of  the  Lagrangian  particle  paths  have  proven  to  be  a  very  useful 
diagnostic  aid  for  determining  the  quality  of  the  solution  and  for  examining  the  physical  mechanism  in 
some  wave  flow  situations.  We  noted  previously  in  the  consideration  of  the  uniform  depth  channel  that 
those  calculations  were  made  with  an  early  version  of  the  SPLISH  code.  While  the  wave  period  and 
wave  drift  were  predicted  accurately  in  that  case,  the  Lagrangian  particle  paths  were  not  as  regular  as 
desired  even  though  the  calculation  was  for  a  simple  channel  flow  without  the  half-cylinder.  Subse¬ 
quent  progress  in  the  development  of  SPLISH  has  substantially  improved  the  quality  of  the  results  for 
particle  paths.  Figure  18  shows  the  Lagrangian  paths  of  three  surface  vertices,  in  Figure  18a  the  vertex 
initially  at  x  —  0  m,  in  18b  the  vertex  at  x  ->  0.6  m,  and  in  18c  the  vertex  at  x  —  1.2m.  The  particle 
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paths  begin  at  step  1 1  after  the  standing  wave  is  fully  formed  and  the  second  pressure  pulse  has  been 
applied.  The  presence  of  the  half-cylinder  in  the  flow  seems  to  affect  the  path  of  the  surface  particle 
above  it  only  slightly.  A  larger  effect  would  be  seen  for  a  larger  wave  amplitude  or  for  vertices  nearer 
to  the  cylinder.  The  particle  paths  in  Figure  18  also  exhibit  a  slight  decrease  in  amplitude  with  time  as 
mentioned  earlier. 

Even  though  there  is  some  loss  of  energy  in  the  numerical  solution  (due  to  incomplete  conver¬ 
gence  of  the  SOR  method,  mentioned  previously),  the  Lagrangian  particle  paths  clearly  show  that,  for 
these  low  wave  reflection  conditions,  SPLISH  is  generating  a  stable,  well-behaved  solution  for  the  sur¬ 
face  wave  flow  over  the  half-cylinder.  The  particle  paths  seem  to  indicate  that  the  wave  reflection  from 
the  half-cylinder  is  sufficiently  low  that  it  neither  adversely  affects  the  quality  of  the  solution  nor 
significantly  changes  the  physical  situation  which  is  being  modeled  numerically.  A  case  in  which  the 
amount  of  wave  reflection  is  significant  and  the  periodic  boundary  conditions  appear  inappropriate  is 
considered  in  the  next  section. 

High  Reflection  Wave  Flow  over  a  Half-Cylinder 

While  the  numerical  results  for  a  low  reflection  case  presented  in  the  preceding  section  are 
promising,  reasonable  numerical  simulations  cannot  be  made  for  high  wave  reflection  from  the  half- 
cylinder  using  periodic  boundary  conditions.  The  results  of  Naftzger  and  Charabarti  (1979)  show  that 
quite  high  reflection  coefficients  can  be  obtained.  Unfortunately,  the  higher  reflection  coefficients  are 
obtained  at  values  of  dl  a  — 1 •  1.0,  i.e.  nearly  complete  blockage  of  the  channel.  With  the  value  of  dl  a 
=  2.0,  as  used  in  the  low  reflection  case,  a  reflection  coefficient,  “  0.24  can  be  obtained  by 
decreasing  ka  (=2ira/\)  from  1.25  to  0.5  (corresponding  to  increasing  \  from  2.5  m  to  6.3  m).  For 
this  value  of  R,  the  reflected  wave  would  have  only  6%  of  the  energy  of  the  incident  wave.  Whereas 
for  d/a  -  1.5,  Rmax  -  0.4  (at  ka  -  0.5)  and  the  reflected  wave  has  16%  of  the  energy  of  the  incident 
wave.  These  latter  values  were  chosen  for  the  high  reflection  case.  The  corresponding  conditions  for 
the  SPLISH  calculations  are  d—  1.0  m,  a  -  0.667  m  and  X  —  8.4  m. 
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In  choosing  these  conditions,  not  only  the  amount  of  an  incident  wave  energy  that  is  reflected  by 
the  bottom  seated  half-cylinder  is  changed,  but  also  the  deep  water/ shallow  water  characteristics  of  the 
wave  flow  were  altered.  For  the  low  reflection  case  the  product  of  the  wave  number  and  depth  has  the 
value  kd  -  2.51  which  gives  almost  deep  water  waves  (kd  >  3.14).  However,  for  the  high  reflection 
case  kd  —  0.75  and  the  wave  flow  is  distinctly  not  in  the  deep-water  regime.  Some  test  calculations  pre¬ 
viously  had  been  made  with  SPLISH  for  a  shallow  water  uniform  depth  channel  with  good  results,  so 
for  the  present  high  reflection  conditions  the  presense  of  the  half-cylinder  was  the  critical  element  in 
the  numerical  simulation. 

Figures  19-22  show  plots  of  the  surface  contour,  the  bottom  contour  and  the  bottom  pressure 
between  x  —  3.2  m  and  x  —  5.2  m  with  the  half-cylinder  centered  at  x  —  4.2m.  The  automatic  grid 
restructuring  algorithms  in  SPLISH  modified  the  initial  grid  slightly  after  the  first  time  step  in  that 
several  triangles  were  added  near  the  intersection  of  the  half-cylinder  and  the  channel  bottom.  Thus  Pb 
for  t  =  0.22  sec  instead  of  for  /“  0.0  sec  is  plotted  in  Figure  19,  closely  approximating  the  initial  still 
water  conditions.  As  in  Figures  8-12  the  pressure  scale  in  kiloPascals  increases  with  depth  and  was 
chosen  so  that  Pb  at  /  =  0  would  correspond  to  the  bottom  contour.  Figure  20,  for  t  —  1.88  sec, 
represents  that  time  in  the  simulation  when  the  water  depth  above  the  half-cylinder  was  at  the  still 
water  level  and  a  wave  crest  was  to  the  left  of  the  half-cylinder.  Figure  21  is  for  t  =  2.68  sec  and  for 
the  wave  crest  directly  above  the  half-cylinder,  and  Figure  22  is  for  t  =  3.28  sec  with  the  wave  crest  to 
the  right  of  the  half-cylinder.  As  before,  the  direction  of  wave  travel  is  from  left  to  right.  These 
figures  show  much  the  same  information  as  did  Figures  8-12,  but  with  the  half-cylinder  larger  and  less 
of  the  surface  wave  visible  in  the  figures  because  of  the  longer  wavelength.  The  curves  for  the  pres¬ 
sure  around  the  half-cylinder  (in  these  figures)  do  not  show  any  detrimental  effect  of  the  wave 
reflection  and  in  fact  appear  to  be  quite  reasonable  representations  of  the  physical  situation. 

The  time  history  of  the  pressure  at  a  point  on  top  of  the  half-cylinder  (9  -  0®)  shown  in  Figure 
23  also  does  not  clearly  indicate  any  detrimental  effect  from  the  reflection  of  the  incident  wave  from 
the  obstacle.  Again  there  is  a  decrease  in  the  amplitude  of  the  numerically  generated  pressure  curve 
much  like  that  seen  in  Figure  17. 
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Somewhat  clearer  evidence  of  the  wave  reflection  is  seen  in  Figures  24  and  25.  Figure  24  shows 
AP/APref  around  the  half-cylinder  at  t  —  2.68  sec  and  t  =  4.14  sec,  corresponding  to  a  crest  and  a 
trough  above  the  half-cylinder.  Figure  25  shows  AP/APREF  at  t  —  1.88  and  3.28  sec  corresponding  to  a 
wave  crest  first  to  the  left  of  the  half-cylinder  and  then  to  the  right.  For  all  four  sets  of  data  the  value 
of  APref  was  the  same  —  the  value  of  A P  at  t  =  2.68  sec  and  at  0  —  0°.  Also  shown  in  these  figures 
are  curves  for  the  fifth-order  theory  (Reference  3)  results  and  curves  for  the  linear  theory  (Reference 
9)  results.  Some  significant  differences  are  noted  between  the  linear  and  fifth  order  theory  results, 
especially  on  the  sides  of  the  half-cylinder.  For  the  times  when  the  crest  or  the  trough  is  above  the 
half-cylinder,  the  fifth  order  theory  results  are  symmetric  with  respect  to  0  —  0°.  That  is,  the  pressure 
fluctuation  distribution  on  the  upstream  (left)  side  of  the  half-cylinder  is  the  same  as  that  on  the  down¬ 
stream  (right)  side.  In  contrast,  for  the  linear  theory  results  there  is  a  pronounced  bulge  in  the  distri¬ 
bution  curve  on  the  upstream  side  of  the  half-cylinder  for  both  crest  and  trough  passage.  Thus,  while 
the  fifth-order  theory  results  would  predict  no  side  force  on  the  half-cylinder  at  either  crest  or  trough 
passage,  the  linear  theory  results  would  show  a  rather  large  force  to  the  right  at  crest  passage  and  a 
comparable  force  to  the  left  at  the  time  that  the  trough  passes  over  the  half-cylinder. 

Another  difference  in  the  results  of  the  linear  and  fifth-order  theories  is  seen  in  the  magnitudes  of 
the  pressure  fluctuation  distributions  at  the  time  of  the  trough  passage.  The  linear  theory  distribution 
curve  has  the  same  magnitude  as  the  curve  for  the  wave  crest  passage,  whereas,  for  the  trough  passage 
over  the  half-cylinder,  the  fifth-order  theory  results  show  a  reduction  in  magnitude  of  the  pressure 
fluctuation  distribution  curve.  For  crest  passage,  the  linear  theory  results  seem  to  more  properly  show 
the  pressure  distribution,  whereas,  for  the  trough  passage,  each  theory  seems  to  model  one  aspect  of 
the  distribution  curve.  The  linear  theory  properly  shows  the  effect  of  wave  scattering  from  the  half¬ 
cylinder,  and  the  fifth-order  theory  properly  shows  the  nonlinear  effect  of  finite-amplitude  waves. 

For  the  pressure  fluctuation  distributions  shown  in  Figure  25  there  is  only  a  slight  difference 
between  the  linear  and  fifth-order  theory  results.  Both  the  vertical  and  horizontal  forces  acting  on  the 
half-cylinder  would  be  slightly  greater  if  calculated  from  the  linear  theory  instead  of  from  the  fifth-order 
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theory.  This  difference  between  the  results  of  the  two  theoretical  approaches  is  quite  small  compared  to 
the  differences  at  the  time  of  crest  and  trough  passage  as  shown  in  Figure  24. 

For  this  case,  the  agreement  is  much  poorer  between  the  theoretical  results  and  the  SPUSH 
results.  For  the  crest  over  the  half-cylinder  (at  t  =  2.68  sec),  the  SPLISH  results  agree  quite  well  with 
the  fifth-order  results,  but  the  linear  theory  results  are  the  more  realistic  because  the  boundary  condi¬ 
tion  on  the  reflected  wave  is  satisfied  in  the  latter  case.  For  the  trough  passage  over  the  half-cylinder 
(at  /  —  4.14  sec)  the  SPLISH  results  seem  to  show  both  an  effect  of  wave  reflection  and  an  effect  of 
finite-amplitude  waves.  However,  a  corresponding  effect  of  wave  reflection  should  have  been  evident 
in  the  distribution  for  the  crest  passage  (at  t  =  2.68  sec)  but  was  not.  Thus  what  appears  to  be  an 
effect  of  wave  reflection  could  instead  be  due  to  the  changing  nature  of  the  wave  flow  which  is  con¬ 
sidered  later  in  the  discussion  of  the  particle  paths  for  this  case. 

The  SPLISH  results  for  t  *=  1.88  and  3.28  sec  shown  in  Figure  25  have  two  interesting  features; 
one  seems  to  show  an  effect  of  wave  reflection,  the  other  almost  certainly  is  due  to  the  effect  of  wave 
reflections  from  neighboring  half-cylinders  in  the  numerical  simulation.  The  SPLISH  results  for  the 
wave  to  the  left  of  the  half-cylinder  (t  =  1.88  sec)  agree  well  with  the  theoretical  results  except  on  the 
downstream  (right)  side  for  0  >  30 °.  This  might  be  an  effect  of  wave  scattering  which  the  theoretical 
results  do  not  show,  but  such  is  not  certain.  The  difference  in  the  theoretical  and  SPLISH  values  of 
AjP/APREf  for  9  >  30°  does  have  the  effect  of  decreasing  the  horizontal  force  on  the  half-cylinder. 
This  effect  is  also  seen  in  Figure  26.  For  the  wave  to  the  right  of  the  half-cylinder  (t  =  3.28  sec)  the 
SPLISH  results  agree  poorly  with  the  theoretical  values.  For  this  time  in  the  wave  cycle,  the  pressure 
distribution  around  the  half-cylinder  from  the  SPISH  results  gives  a  value  for  the  horizontal  force 
which  is  much  too  small  as  is  discussed  further  later.  This  poor  agreement  between  the  numerical  and 
theoretical  results  seems  to  be  caused  by  the  use  of  the  periodic  boundary  conditions  in  the  SPLISH 
code  with  the  resulting,  non-physical,  wave  reflections  from  the  neighboring  half-cylinders  which  exist 
only  in  the  numerical  simulation.  The  interaction  between  the  periodic  boundary  conditions  and  the 
wave  reflections  is  discussed  further  in  the  consideration  of  the  Lagrangian  particle  paths. 
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Time-history  plots  of  the  numerically-calculated  force  components  on  the  half-cylinder  are  shown 
in  Figures  26  and  27.  Again,  as  in  Figures  15  and  16,  the  plotting  of  the  curves  is  begun  after  the  ini¬ 
tial  transients  die  out.  Linear  wave  calculations  from  Naftzger  and  Chakrabarti’s  predictions  (1979)  are 
shown.  As  for  the  previous  low  reflection  case,  the  vertical  force  component  curve  shown  in  Figure  27 
appears  quite  similar  to  the  pressure  curve  for  the  location  on  top  of  the  half-cylinder  shown  in  Figure 
23.  The  local  extrema  of  the  two  curves  show  quite  similar  characteristics.  Indeed,  both  the  force  and 
the  pressure-time  variations  seem  most  nearly  like  slightly  modified  sine  or  cosine  curves. 

Such  is  not  the  case,  however,  for  the  horizontal  force  component  Fx  as  shown  in  Figure  26. 
There  is  part  of  a  sine  curve  from  z  =  1  to  2  sec,  then  something  like  a  saw-tooth  curve  from  t  =  2  to 
5  sec,  followed  by  a  full  cycle  of  a  sine  wave  and  then  another  saw-tooth.  Consideration  was  given  to 
the  possibility  that  this  behavior  was  due  to  the  difficulty  in  accurately  calculating  the  horizontal  force 
component  along  the  steep  vertical  sides  of  the  half-cylinder.  However,  this  behavior  of  the  Fx  curve 
did  not  appear  in  Figure  16  for  the  low  reflection  case,  and  so  the  manner  of  calculating  Fx  did  not 
appear  to  be  the  cause  of  the  irregularities  in  the  horizontal  force  component.  Considering  Figures  25 
and  26  together  shows  that  at  t  =  1.88  seconds  the  pressure  distribution  from  the  numerical  simulation 
agreed  well  with  the  theoretical  distribution  and  that  the  horizontal  force  component  still  is  quite  regu¬ 
lar.  However,  at  t  =  3.28  sec,  the  numerical  and  theoretical  pressure  distributions  agree  poorly.  Closer 
examination  of  the  two  curves  in  Figure  25  shows  that  they  give  similar  values  for  Fy  but  very  different 
values  for  Fx.  Indeed,  as  may  be  seen  in  Figure  26,  the  magnitude  of  Fx  at  t  —  3.28  sec  is  much  lower 
than  would  be  obtained  if  Fx  followed  the  expected  sine  curve.  This  behavior  of  the  horizontal  force  — 
time  history  can  be  attributed  to  the  interaction  between  the  reflected  waves  and  the  periodic  boundary 
conditions.  Further  effects  of  this  interaction  are  seen  in  the  Lagrangian  particle  paths  which  are  con¬ 
sidered  next. 

In  Figures  28-31  the  Lagrangian  particle  paths  are  shown  for  surface  vertices  initially  at  x  =  0.0 
m,  2.0  m,  4.2  m  and  6.4  m.  The  initial  positions  of  the  particles  on  the  still  water  level  are  indicated  by 
a  circled  x.  The  anomalies  in  the  Lagrangian  particle  paths  are  the  most  striking  features  of  these 
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figures.  In  addition  to  specific  comments  about  the  individual  figures,  several  general  observations  can 
be  made.  Since  for  these  conditions  the  progressive  waves  can  be  characterized  as  approaching  shallow 
water  waves,  the  particle  paths  should  be  more  elliptic  in  nature  than  circular.  This  is  the  case  in  Fig¬ 
ures  29  and  31  for  the  vertices  at  x  =  2.0  m  and  6.4  m.  However,  in  Figures  28  and  30  the  particle 
paths  for  the  vertices  at  x  =  0.0  m  and  4.2  m  are  much  more  circular  in  nature,  at  least  initially.  The 
overall  effect  is  that  at  x  =  0.0  m  and  4.2  m  the  wave  amplitude  is  approximately  twice  what  it  is  at  x 
“  2.0  and  6.4  m. 

In  all  four  figures  it  appears  that  the  open  orbits  of  the  particle  paths  are  collapsing  to  flattened, 
inclined,  paths  along  the  streamlines  characteristic  of  a  standing  wave  pattern.  It  is  here  that  the 
interactions  of  the  wave  reflections  and  periodic  boundary  conditions  are  most  evident.  In  the  numeri¬ 
cal  simulation  the  bottom  seated  half-cylinder  seen  in  the  figures  represents  one  member  of  an  array  of 
half-cylinders,  uniformly  spaced  one  wavelength  apart.  The  periodic  boundary  conditions  require  that 
conditions  at  x  —  0  m  and  at  x  =  8.4  m  be  the  same,  so  that  the  computational  domain  is  one  of  a  set 
of  identical  domains.  Because  of  the  high  reflection  coefficients  in  this  physical  situation,  the  reflected 
waves  from  one  barrier  should  interfere  with  the  transmitted  waves  from  the  preceding  barrier.  If  a 
progressive  wave  is  viewed  as  a  superposition  of  standing  waves,  interference  would  be  expected  among 
the  standing  wave  components  of  the  reflected  and  transmitted  waves.  Thus,  with  repeated  reflections 
of  waves  between  adjacent  half-cylinders,  the  initial  progressive  wave  seems  to  transform  into  standing 
waves.  Also,  the  manner  in  which  the  progressive  wave  is  initiated  does  not  represent  an  incoming 
wave  field,  but  rather  a  wave  field  suddenly  springing  up  (under,  of  course,  the  influence  of  the  applied 
pressure  pulses). 

We  note  that  at  x  -  4.2  m  the  surface  vertex  is  a  full  wave  length  away  from  the  neighboring 
half-cylinders  (excepting  the  one  directly  under  it),  and  it  should  move  independently  of  the  effects  of 
reflected  waves  for  a  full  wave  cycle  or  more.  This  does  indeed  seem  to  be  the  case,  as  shown  in  Fig¬ 
ure  30.  In  a  similar  manner,  the  particle  at  x  —  0  m  should  be  free  of  the  effect  of  wave  reflection  for 
three  quarters  of  a  wave  cycle  as  the  trough  formed  at  x  -  2.1  m  reaches  the  half-cylinder  at  774  sec 
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and  the  first  reflection  reaches  x  =  0  in  a  time  of  772  sec  or  at  t  —  3774  sec.  As  shown  in  Figure  28, 
this  seems  to  be  the  case.  The  particle  at  x  —  2.1  m  should  be  first  influenced  by  the  wave  reflection  at 
772  sec  as  the  trough  takes  7/4  sec  to  reach  the  half-cylinder  and  the  reflection  takes  7/4  sec  to  return. 
The  particle  path  in  Figure  29  suggests  this  sort  of  behavior. 

The  Lagrangian  particle  path  for  the  vertex  initially  at  x  =  6.4  m  should  show  about  one  wave 
cycle  of  motion  before  being  affected  by  any  wave  reflection.  The  trough  initially  formed  at  x  •=  10.5 
m  (i.e.,  8.4  m  +  A/4)  should  reach  the  half-cylinder  to  the  right  in  7/4  sec  and  the  reflection  back 
should  reach  x  —  6.4  m  in  an  additional  37/4  sec.  However,  there  is  some  effect  on  the  particle  path  at 
7/2  sec,  as  observed  in  the  bent  shape  of  the  ellipse,  and  then  nearly  a  full  wave  cycle  passes  before 
the  particle  path  appears  to  collapse  to  that  of  a  standing  wave.  The  first  irregularity  could  be  due  to 
the  half-cylinder  having  some  unrecognized  effect  on  the  passage  of  the  first  wave  trough  over  it. 

Most  of  the  characteristics  of  the  Lagrangian  particle  paths  in  Figures  28-31  seem  to  be  due  to  the 
reflections  of  waves  from  adjacent  half-cylinders.  However,  one  characteristic  remains  unexplained. 
The  wave  heights  at  x  =  A/4  and  3A/4  are  only  half  of  the  wave  heights  at  x  -  0  and  A/2.  Wave  rein¬ 
forcement  and  cancellation  does  not  provide  a  fully  satisfactory  explanation.  The  points  at  x  —  0  and 
A/2  are  points  where  the  reflected  wave  should  reinforce  the  incident  wave.  Likewise,  the  points  at 
x  =  A/ 4  and  3A/ 4  are  points  where  the  reflected  waves  should  partially  cancel  the  incident  waves.  As 
the  initial  wave  is  formed  by  the  first  pressure  pulse,  the  nodes  are  at  x  —  0  m  and  4.2  m  (i.e.,  at  0  and 
A/2),  the  trough  is  at  x  =  2.1  m  (A/4),  and  the  crest  is  at  x  =*  6.3  m  (3A/4).  However,  at  least  half  of 
a  wave  cycle  (if  not  full  wave  cycle)  should  pass  before  any  such  wave  reinforcement  or  cancellation 
occurs.  The  particle  paths  in  Figures  28  and  30  show  that  the  wave  reaches  full  amplitude  at  7/4  sec 
(measured  from  the  second  pressure  pulse  which  forms  the  progressive  wave)  with  no  assistance  of 
reinforcement  from  a  reflected  wave.  The  particle  paths  in  Figures  28-31  almost  seem  to  indicate  that 
the  primary  wave  flow  is  modulated  by  a  secondary  standing  wave.  Unless  the  half-cylinder  poses  more 
blockage  to  the  progressive  wave  flow  than  had  been  assumed,  it  is  not  clear  what  the  source  of  such  a 
standing  wave  component  might  be. 
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While  the  development  of  a  standing  wave  is  perhaps  most  obvious  in  Figures  28-31,  as  discussed 
above,  other  figures  also  suggest  such  an  occurrence.  For  example,  the  pressure-time  history  for  the 
top  of  the  half-cylinder  (shown  in  Figure  23)  is  really  quite  remarkably  regular  and  in  good  agreement 
with  fifth-order  wave  theory  considering  the  irregularities  in  the  flow  field.  Yet  this  pressure-time  his¬ 
tory  would  be  quite  consistent  with  a  standing  wave  flow  field  with  crests  and  troughs  nearly  above  the 
half-cylinder.  The  vertical  force  time  history  shown  in  Figure  27  is  also  remarkably  regular  but  is  also 
consistent  with  the  existence  of  a  standing  wave.  The  irregularities  in  the  horizontal  force  shown  in 
Figure  26  could  be  explained  in  part  by  a  standing  wave.  That  is,  with  crests  and  troughs  of  a  standing 
wave  centered  above  the  half-cylinder,  there  would  be  little,  if  any,  horizontal  force  on  the  half¬ 
cylinder.  Indeed,  the  horizontal  force  calculated  in  the  numerical  simulation  is  much  less  than 
predicted  by  the  linear  wave  scattering  theory  of  Naftzger  and  Chakrabarti  (1979). 

Also,  pressure  distributions  around  the  half-cylinder,  such  as  those  shown  in  Figures  24  and  25, 
for  later  times  in  the  SPLISH  calculation  show  a  further  effect  of  a  developing  standing  wave. 
Specifically,  for  times  corresponding  to  crests  and  troughs,  the  distribution  of  A P/APREF  are  much  like 
those  shown  in  Figure  24.  However,  the  distribution  for  intermediate  times  (as  in  Figure  25)  show  an 
increasingly  large  difference  between  the  SPLISH  results  and  the  results  from  wave  theory.  In  particu¬ 
lar,  the  SPLISH  results  approach  AP/APREF  =  0  even  more  strongly  than  was  the  case  for  t  —  3.28  sec. 
Thus,  at  times  when  the  wave  theories  would  give  maximum  side  forces  on  the  half-cylinder,  the 
SPLISH  results  show  only  very  small  side  forces.  Such  behavior  is  fully  consistent  with  a  standing 
wave  pattern  with  crests  and  troughs  approximately  above  the  half-cylinder.  While  we  see  considerable 
evidence  for  a  standing  wave  pattern  in  the  numerical  results,  most  of  the  evidence  points  to  the  gra¬ 
dual  development  of  the  standing  wave  after  repeated  wave  reflections.  It  remains  unexplained  why  in 
the  SPLISH  results  the  wave  height  at  x  -  0  m  and  at  x  -  4.2  m  is  nearly  twice  the  wave  height  at  x  — 
2.0  m  and  x  -  6.4  m. 

The  results  for  the  low  reflection  case  indicate  that,  even  with  the  periodic  boundary  conditions, 
SPLISH  can  be  used  with  reasonable  confidence  for  calculating  the  flow  of  waves  over  obstacles  (thus 
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wave-structure  interactions)  if  wave  reflections  are  not  significant.  However,  for  cases  with  significant 
wave  reflection  in  the  flow,  alternate  (non-periodic)  boundary  conditions  are  clearly  needed,  since 
under  the  periodic  boundary  conditions  a  different  physical  situation  is  modeled  with  a  correspondingly 
different  flow  as  a  result.  Further  alternate  techniques  for  initiating  the  progressive  wave  seem  to  be 
needed.  It  would  also  be  very  desirable  to  achieve  large  reductions  in  computing  costs  and  time.  Each 
of  the  SPLISH  calculations  for  the  two  sets  of  conditions  presented  here  required  about  20  minutes  of 
computing  time  on  the  Texas  Instruments  (TI)  Advanced  Scientific  Computer  (ASC)  at  NRL  at  a  cost 
of  approximately  $400.  These  times  and  costs  were  incurred  even  though  the  code  had  been  especially 
developed  to  take  maximum  advantage  of  the  vector  characteristics  of  the  ASC.  Thus,  significant 
reductions  in  computing  costs  and  times  may  be  difficult  to  achieve.  Additionally,  the  large  memory 
requirements  of  the  code  would  make  it  difficult  to  implement  SPLISH  on  other  computers. 

In  the  above  discussion  we  have  only  compared  the  results  obtained  with  SPLISH  with  results 
from  classical  linear  and  fifth-order  wave  theories.  Even  though  most  of  these  comparisons  showed 
good  to  very  good  agreement,  it  is  also  appropriate  to  compare  all  of  the  calculations  with  experimental 
data.  This  is  done  in  the  next  section. 

Comparisons  With  Experimental  Data 

Several  wave  flow  experiments  were  performed  in  the  NRL  wave  channel  with  a  bottom-seated 
half-cylinder  so  that  actual  pressure  measurements  could  be  compared  both  to  the  numerical  results  and 
to  the  models  developed  by  Chakrabarti  and  Naftzger  (1974,  1979).  A  1.07  m  (3-1/2  ft)  diameter 
cylinder,  which  spanned  the  entire  width  of  the  channel,  was  placed  about  one-half  of  the  channel’s 
length  from  the  mechanical  wavemaker.  At  the  other  end  of  the  wave  channel  a  sloping,  porous  beach 
with  a  rubberized  horsehair  blanket  served  to  absorb  nearly  all  of  the  incident  wave  energy.  A  more 
detailed  description  of  the  facility  is  available  (see  reference  Cl  of  Appendix  C).  Nineteen  equally- 
spaced  (A 9  —  10°)  pressure  taps  were  located  around  the  circumference  of  the  bottom-seated  half¬ 
cylinder  at  its  midsection.  The  individual  taps  were  connected  to  a  differential  pressure  transducer  by  a 
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rotary  pressure  switch.  The  wave  height  along  the  channel  was  obtained  from  several  traversing 
capacitance-type  wave  gauges.  Calibrations  were  performed  on  all  sensors  before  and  after  each  test 
series  to  insure  that  the  overall  accuracy  of  the  measurement  system  remained  well  within  ±5  percent. 
The  pressure  and  waveheight  signals  were  digitized  and  processed  by  means  of  a  Hewlett-Packard 
5420A  Digital  Signal  Analyzer.  The  experimental  systems  and  methods  are  discussed  further  in  Appen¬ 
dix  C. 

The  cylinder  reflection  coefficient  R,  defined  as  the  ratio  of  reflected  to  incident  wave  amplitudes, 
was  obtained  for  a  range  of  wavelengths  and  water  depths.  The  results  of  these  measurements,  after  a 
rudimentary  adjustment  for  the  unwanted  effects  discussed  in  Appendix  C,  are  plotted  in  Figure  32 
versus  the  wavenumber  ka  =  2 ira/X.  and  for  several  values  of  the  relative  water  depth  d/a.  Also 
shown  in  the  figure  (as  solid  lines)  are  theoretical  values  based  on  linear  theory  with  wave  scattering 
which  were  taken  from  the  paper  by  Naftzger  and  Chakrabarti  (1979)  and  are  in  general  agreement  with 
the  observed  cylinder  reflections.  It  should  be  noted,  however,  that  some  care  was  taken  to  avoid  finite 
amplitude  effects  through  the  use  of  small  wave  steepness  ratios  (///A  <  0.05).  The  wave  steepness 
values  during  the  d/a  —  1.25  tests  were  reduced  further,  typically  less  than  0.02,  in  order  to  avoid 
second  harmonic  wave  generation  at  the  cylinder.  This  nonlinear  effect  seemed  to  be  associated  with 
the  finite  waveheight  being  a  significant  fraction  of  the  finite  water  depth  at  the  top  of  the  cylinder. 
The  second  harmonic  wave  could  be  seen  as  the  formation  of  a  secondary  crest  when  a  trough  was  over 
the  cylinder.  A  rule  of  thumb  emerged  from  this  observation  and  from  several  wave  gauge  spectral 
records  (FFT)  wherein  second  harmonic  amplification  did  not  occur  for  waveheights  less  than  about 
1/7  to  1/10  of  the  water  depth  over  the  cylinder.  This  criterion  for  the  validity  of  the  analytical  method 
is  considerably  more  stringent  than  the  one  originally  proposed  by  Naftzger  and  Chakrabarti.  The  onset 
and  form  of  this  nonlinearity  will  be  the  subject  of  further  investigation  in  the  near  future. 

From  the  range  of  water  depths  and  wavelengths  shown  in  Figure  32,  two  cases  were  selected  for 
detailed  pressure  studies  and  for  comparisons  with  the  theoretical  and  numerical  results.  The  first  is  a 
relatively  low  reflection  case  {d/a  —  2.0  and  ka  ■■  1.25),  so  that  R  <  0.05.  This  represents  a  situation 
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where  good  agreement  between  the  three  sets  of  results  was  expected.  The  second  case  is  a  relatively 
high  reflection  condition  (d/a  =■=  1.5,  ka  =  0.5,  R  =  0.4)  which  was  selected  as  a  significant  test  for 
SPLISH  with  regard  to  the  use  of  periodic  boundary  conditions  in  the  code. 

Low  Reflection  Case 

A  comparison  between  the  linear  asymptotic  theory  and  the  experimental  results  is  given  in  Fig¬ 
ure  33.  The  asymptotic  form  was  selected  for  comparison  here  because  of  its  simplicity.  The  ampli¬ 
tudes  and  the  relative  phases  of  the  pressure  fluctuations  are  plotted  against  azimuthal  position  around 
the  cylinder  as  measured  in  degrees  from  the  top  of  the  cylinder.  All  pressures  are  nondimensionalized 
by  the  maximum  pressure  amplitude  predicted  by  the  theory.  This  maximum  occurs  at  the  top  of  the 
cylinder  in  this  case  (9  ■*  0T  The  phase  angles  are  also  plotted  relative  to  the  linear  theory  value  at 
0  =  0°.  The  two  sets  of  experimental  data  correspond  to  two  wave  steepnesses  and  indicate  the  absence 
of  finite  amplitude  effects  in  this  range.  Except  for  the  data  at  the  largest  azimuthal  angles,  the  agree¬ 
ment  between  the  various  results  is  reasonably  good  and  supports  the  use  of  the  linear  asymptotic 
theory.  The  differences  which  exist  at  the  azimuthal  extremes  are  possibly  a  residual  of  the  reflection 
adjustment  and/or  a  consequence  of  a  small  gap  which  existed  between  the  bottom  of  the  cylinder  and 
the  wave  tank  floor. 

Since  the  numerical  results  are  Lagrangian  whereas  the  theoretical  and  experimental  results  are 
Eulerian,  the  simplest  format  for  comparing  all  three  is  the  pressure  distribution  around  the  cylinder  at 
selected  times  in  the  wave  cycle.  Figures  34  and  35  presents  several  such  comparisons  for  the  experi¬ 
mental  conditions  cited  above.  Again  the  three  sets  of  results  compare  well  except  near  the  bottom  of 
the  cylinder.  Although  the  discrepancy  appears  to  be  small  in  such  a  plot  the  effect  on  a  measured  hor¬ 
izontal  component  of  the  wave  force  can  be  significant. 

High  Reflection  Case 

Similar  results  and  comparisons  for  a  relatively  high  reflection  case  (d/a  -  1.5,  ka  -  0.5, 
R  =  0.4)  are  shown  in  Figures  36  through  38.  As  before,  the  pressures  are  normalized  by  the  magni¬ 
tude  at  the  top  of  the  cylinder.  In  this  case,  however,  the  maximum  pressure  fluctuation  occurs  on  the 
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upstream  side  of  the  cylinder  where  both  the  incident  and  reflected  waves  are  present.  The  experimen¬ 
tal  pressures  and  phases  shown  in  the  figures  were  obtained  from  the  data  after  a  correction  for  the 
additional  standing  waves  in  the  laboratory  channel  (see  Appendix  C).  This  correction,  often  large, 
neglected  secondary  reflections  and  was  based  on  linear  wave  theory  which  may  be  responsible  for  all  or 
part  of  the  discrepancies  between  results.  The  theoretical  results  shown  for  the  method  of  Naftzger  and 
Chakrabarti  were  kindly  provided  by  them  (R.A.  Naftzger,  private  communication,  1980). 

In  spite  of  the  difficulties  implied  above  and  discussed  in  Appendix  C,  the  agreement  between  the 
linear  theory  and  the  experiments  is  fair.  The  previously  mentioned  restriction  on  wave  amplitude 
ought  to  be  reiterated.  The  comparisons  with  the  SPLISH  results  show  the  detrimental  effect  of 
surprisingly  good  agreement  considering  the  use  of  periodic  boundary  conditions  in  the  code. 
SUMMARY  AND  CONCLUSIONS 

A  finite-difference  numerical  method  for  solving  the  governing  equations  of  motion  for  inviscid, 
irrotational  flow  with  a  free  surface  using  a  Lagrangian  triangular  grid  has  been  used  and  shown  to  yield 
reasonable  results.  Calculations  for  progressive  surface  wave  flows  have  given  results  for  the  wave 
period,  the  drift  velocity  and  the  surface  particle  movements  with  are  in  good  agreement  with  results 
obtained  from  classical  wave  theory. 

Calculations  for  the  passage  of  waves  over  a  submerged  obstacle  are  encouraging  and  show  some 
promise  of  providing  practical  results  for  more  complex  flows  such  as  breaking-wave  forces.  These  cal¬ 
culations  demonstrate  the  adaptability  which  the  triangular  grid  provides.  The  advantages  of  the 
Lagrangian  formulation  are  shown  in  that  the  grid  conforms  to  the  fluid  area  and  that  no  interpolation 
is  needed  to  locate  the  free  surface  or  the  surface  of  the  submerged  obstacle. 

Two  cases  of  flow  over  a  submerged,  half-cylindrical  obstacle  have  been  considered.  The  results 
for  the  low  wave  reflection  case  indicate  that,  even  with  the  periodic  boundary  conditions,  SPLISH  can 
be  employed  with  reasonable  confidence  to  calculate  the  flow  of  waves  over  obstacles  (thus  wave- 
structure  interactions).  However,  for  cases  with  significant  wave  reflection  in  the  flow,  continuative  or 
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radiative  boundary  conditions  are  clearly  desirable,  if  not  necessary,  since  the  periodic  boundary  condi¬ 
tions  lead  to  the  simulation  of  a  flow  quite  different  from  that  desired.  Also,  alternative  techniques  for 
initiating  the  traveling  surface  wave  seem  to  be  required. 

While  additional  work  on  non-periodic  boundary  conditions,  initial  conditions  and  gridding  of  the 
flow  field  is  clearly  necessary,  the  present  numerical  method  is  a  promising  computational  tool  for 
predicting  the  velocities  and  pressure  fields  produced  by  nonlinear,  finite-amplitude  regular  wave 
motions.  It  is  hoped  that  future  developments  will  extend  the  capabilities  of  codes  such  as  SPLISH  to 
transient  phenomena  such  as  wave  breaking  and  wave  slamming  and  that  significant  reductions  in  com¬ 
puting  times  and  costs  can  be  achieved.  Reduction  in  the  memory  requirements  would  also  be  desir¬ 
able. 
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Figure  1  —  An  initial  triangular  mesh 
for  the  SPL1SH  computer  code. 
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Figure  2  —  An  example  of  a  SPLISH-generated 
progressive  wave. 
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Figure  3  —  An  example  of  a  SPLISH-generated  computational 
grid  with  triangle  reconnections. 
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Figure  4  —  A  comparison  between  Stokes’  second-order  wave  theory 
and  SPLISH-generated  particle  paths  for  wavelength  X  —  2.5  m  and 
depth  d  -  1  m. 


Figure  6  —  An  example  of  a  SPLISH-generated  triangular  grid  for 
computing  the  wave  flow  over  a  bottom  seated  half-cylinder. 
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half-cylinder  at  t  -  0  sec  (time  step  —  1,  wavelength  V  -  2.5  m,  depth  d  -  1  m,  radius  a  -  0.5  m). 
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Figure  10  —  The  surface  contour  and  the  bottom  pressure  from  a  numerical  calculation  of  the  wave  flow  over  a  bottom  seated 
half-cylinder  at  t  -  1.28  sec  (time  step  -  33,  wavelength  X  ~  2.5  m,  depth  d  —  1  m,  radius  a  0.5  m). 


half-cylinder  at  t  -  1.64  sec  (time  step  -  43,  wavelength  X  -  2.5  m,  depth  d  m,  radius  a  0.5  m). 
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Figure  14  —  A  comparison  of  the  distribution  of  pressure  fluctuations  around  the  half-cylinder  from  fifth-order  wave  theory  and 
from  a  SPLISH  calculation,  for  wave  crest  to  the  left  and  the  right  of  the  half-cylinder  (wavelength  \  -  2.5  m,  depth  d  -  1  m, 
radius  a  —  0.5  m). 
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component  of  force  on  the  half-cylinder  from  a  SPLISH  calculation  (wavelength  X  -  2.5  m,  depth 
d  -  1  m,  radius  a  -  0.5  m). 
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Figure  17  —  A  comparison  of  the  pressure  time  history  at  the  top  of  a  bottom  seated  half-cylinder  from  a  SPLISH  calculation 
and  from  fifth-order  wave  theory  (wavelength  A  —  2.5  m,  depth  d  «  1  m,  radius  a  ™  0.5  m,  wave  amplitude  h  ™  0.38  m). 
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Figure  18  —  Particle  paths  on  the  free  surface  from  a  numerical  calculation  of  wave  flow  over  a  bottom 
seated  half-cylinder  (wavelength  A  -  2.5  m,  depth  d  —  1  m,  radius  a  —  0.5  m). 
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Figure  19  —  The  surface  contour  and  the  bottom  pressure  from  a  numerical  calculation  of.  the  high  reflection  wave  flow  over  a 
bottom  seated  half-cylinder  at  t  -  0.2  sec  (time  step  -  11,  wavelength  X  -  8.4  m,  depth  d  -  1  m,  radius  a  0.667  m). 
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e  surface  contour  and  the  bottom  pressure  from  a  numerical  calculation  of  the  high  reflection  wave  flow  over  a 
half-cylinder  at  t  -  2.68  sec  (time  step  -  135,  wavelength  X  -  8.4  m,  depth  d  —  1  m,  radius  a  —  0.667  m). 
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comparison  of  the  pressure-time  history  at  the  top  of  a  bottom  seated  half-cylinder  from  a  SPLISH  calculation 
rave  theory  calculation  for  high  reflection  wave  flow  (wave  length  \  -  8.4  m.  depth  d-  1  m.  radius  a  -  0.667  m 


THEORY;  (REFERENCE  9  -  LINEAR;  REFERENCE  3  -  FIFTH  ORDER) 

-  REF  9  -  CREST  OVER  o  SPLISH  -  CREST  OVER 

- REF  3  -  CREST  OVER  T  -  2.68  SEC 

-  REF  9  -  TROUGH  OVER  □  SPLISH  -  TROUGH  OVER 

-  REF  3  -  TROUGH  OVER  T  -  4.14  SEC 
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Figure  24  --  A  comparison  of  the  distribution  of  pressure  fluctuations  around  the  half-cylinder  from  linear  (Reference  9)  and 
fifth-order  (Reference  3)  wave  theory  and  a  SPLISH  calculation,  for  wave  crest  and  trough  passage  over  the  cylinder  (high 
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(ure  25  —  A  comparison  of  the  distribution  of  pressure  fluctuations  around  the  half-cylinder  from  linear  (Reft 
"th-order  (Reference  3)  wave  theory  and  a  SPLISH  calculation  for  wave  crest  to  the  left  and  to  the  right  of  the 
iigh  reflection  wave  flow;  wavelength  \  -  8.4  m,  depth  d  -  1  m,  radius  a  -  0.667  m). 
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Figure  2$  -  The  particle  path  (at  x  -  0  m)  on  the  free  surface  from  a  SPLISH  calculation  of  the  high  reflection  wave  flow  over  a 
bottom  seated  half-cylinder  (wavelength  \  —  8.4  m,  depth  d  —  1  m,  radius  a  “  0.667  m). 


—  2.0  m)  on  the  free  surface  from  a  SPLISH  calculation  of  the  high  reflection  wave 
half-cylinder  (wavelength  \  ~  8.4  m,  depth  d  ■*  1  m,  radius  a  ~  0.667  m). 


Figure  31  —  The  particle  path  (at  x  —  6.4  m)  on  the  free  surface  from  a  SPLISH  calculation  of  the  high  reflection  wave  flow 
flow  over  a  bottom  seated  half-cylinder  (wavelength  X  -  8.4  m,  depth  <f  —  1  m,  radius  a  -  0.667  m). 
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Figure  32.  Experimental  (•)  and  theoretical  (— )  values  for  the  wave  reflection  coefficient  of  a  bottom-seated  half  cylinder. 
Theoretical  values  are  from  Naftzger  and  Chakrabarti  (1979). 
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Figure  33.  A  comparison  of  theoretical  (Chakrabarti  and  Naftzger,  1974)  and  experimental  results  for  the  azimuthal  distribution 
of  pressure  fluctuation  magnitudes  and  relative  phases  on  a  bottom-seated  half  cylinder  in  waves  (low  reflection  wave  flow,  d! a 


EXP  -  CREST  OVER 
SPLISH  -  CREST  OVER 


EXP  -  WAVE  TO  LEFT 
SPLISH  -  HAVE  TO  LEFT 


re  35.  A  comparison  of  theoretical  (Chakrabarti  and  Naftzger,  1974),  numerical,  and  experimental  results  for  selected  ii 
aneous  pressure  distributions  on  a  bottom-seated  half-cylinder  in  waves,  for  wave  crest  to  the  left  and  the  right  of  the  hal 
iter  (low  reflection  wave  flow,  dt  a  —  2.0,  ka  —  1.25,  R  <  0.05). 
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Figure  36.  A  comparison  of  theoretical  (Naftzger  and  Chakrabarti,  1979)  and  experimental  results  for  the  azimuthal  distribution 
of  pressure  fluctuation  magnitudes  and  relative  phases  on  a  bottom-seated  half  cylinder  in  waves  (high  reflection  wave  flow,  d/a 
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Appendix  A 

NUMERICAL  METHODS  AND  ALGORITHMS  FOR  TRIANGULAR,  LAGRANGIAN  GRIDS 
Numerical  Methods 

A  Lagrangian  treatment  is  particularly  attractive  due  to  the  nature  of  the  fluid  flows  of  interest. 
For  wave-body  interaction  problems  the  free  surface  may  dominate  the  evolution  of  the  flow,  and  the 
submerged  or  partially  submerged  bodies  may  move  as  well,  so  that  a  simple  and  accurate  treatment  of 
the  flow  near  boundaries  is  essential.  The  strength  of  the  Lagrangian  formulation  resides  in  the  fact 
that  fluid  elements  are  advected  with  the  flow.  Therefore  the  grid  points  which  define  surfaces  remain 
at  those  surfaces  and  permit  the  maximum  accuracy  in  formulation  of  boundary  conditions.  The  varia¬ 
tion  of  grid  resolution  which  must  occur  in  Eulerian  schemes  as  surfaces  pass  through  cells  is  also 
alleviated  and  the  numerical  diffusion  across  boundaries  minimized.  Most  importantly  the  nonlinear 
convective  terms  are  not  present  in  the  Lagrangian  formulation,  resulting  in  higher  accuracy  and  less 
stringent  resolution  requirements. 

However,  the  great  strength  of  the  Lagrangian  approach  is  also  its  greatest  weakness.  The  advec- 
tion  of  the  mesh  with  the  flow  leads  to  large  mesh  deformations  and  a  corresponding  decrease  in  accu¬ 
racy  in  both  finite  difference  and  finite  element  methods.  We  will  first  discuss  the  resolution  of  the 
problems  associated  with  deformed  grid  resolution  in  standard  Lagrangian  computer  codes.  The 
deficiencies  in  earlier  solutions  are  examined  and  an  alternative  solution  is  proposed  and  described.  A 
detailed  description  of  its  properties,  implementation  and  tests  of  accuracy  are  described  in  the  next 
section. 

Deformed  Lagrangian  Grids 

To  illustrate  the  effect  of  grid  deformation  let  us  examine  first-  or  second-order  accurate  finite- 
difference  methods.  The  mesh  points  commonly  used  to  evaluate  gradients  and  Laplacians  are  shown 
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in  Figure  Ala  for  a  regular  grid.  Figure  Alb  illustrates  a  simple  and  common  grid  distortion  produced 
by  shear  flow.  A  well-formulated  Lagrangian  finite-difference  algorithm  will  properly  account  for  the 
angles  between  grid  lines  and  the  variable  mesh  spacing  produced  by  this  distortion.  Nevertheless, 
numerical  approximations  based  on  this  mesh  can  still  be  grossly  in  error  because  differences  no  longer 
involve  neighboring  vertices.  Mesh  points  now  closer  to  the  central  vertex  do  not  enter  into  the 
approximation,  while  those  further  removed  do.  As  shown  in  Figure  A2,  higher  order  approximations 
may  lead  to  even  greater  error.  Figure  A2a  shows  the  vertices  commonly  used  in  higher-order  approxi¬ 
mations.  Figure  A2b  illustrates  that  these  approximations  on  a  distorted  mesh  may  include  vertices 
which  are  even  further  removed  from  the  central  vertex  while  neglecting  other  vertices  which  lie  closer. 
In  other  words  the  distorted  mesh  cannot  be  used  to  self-consistently  improve  the  accuracy  of  the 
approximation.  The  problem  can  be  resolved  only  by  differencing  over  the  appropriate  vertices.  That 
is,  the  mesh  distortion  must  be  reduced.  This  same  conclusion  holds  true  for  the  finite  element 
method.  Whether  triangular  or  quadrilateral  elements  are  chosen  the  accuracy  on  the  deformed  grid 
must  be  diminished  regardless  of  the  basis  functions,  simply  because  the  grid  joins  the  wrong  vertices. 

Rezoning 

The  traditional  solution  to  this  problem  has  been  to  perform  an  Eulerian  rezoning  phase  which 
allows  the  mesh  to  pass  through  the  fluid.  For  example,  on  a  quadrilateral  mesh  a  central  vertex  may 
be  moved  to  the  average  position  of  its  neighbors.  The  Eulerian  motion  of  the  vertex  causes  fluid  to 
pass  through  the  four  quadrilateral  cell  sides  meeting  at  that  vertex,  affecting  the  physical  quantities 
specified  in  all  four  cells  and  the  four  neighboring  vertices  as  well  as  the  central  vertex.  The  result  of 
the  rezone  phase  is  therefore  the  introduction  of  artificial  diffusion  over  one  grid  spacing  in  any  direc¬ 
tion.  A  related  but  separate  question  arises  over  the  accuracy  of  the  rezone  algorithm.  The  rezone 
phase  constitutes  a  separate  Eulerian  fluid  flow  calculation.  In  fluid  flows  with  shear  or  near  stagnation 
points  the  vertex  motion  accomplished  during  the  rezone  phase  can  be  equal  (and  opposite)  to  the 
advective,  Lagrangian  vertex  motion.  In  such  cases  the  flow  calculation  is  Eulerian.  Unless  the  algo¬ 
rithms  adopted  for  rezoning  approach  the  accuracy  available  for  finite-difference  approximations  of  the 
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advective  terms  on  an  Eulerian  mesh,  the  accuracy  of  the  calculation  is  necessarily  degraded  and  even 
more  diffusive  than  a  purely  Eulerian  calculation.  A  third  objection  to  rezoning  is  that  it  is  useless  in 
fluid  flows  which  force  changes  from  simply-connected  to  multiply-connected  regions,  as  in  the  case  of 
wave-breaking  or  bubble  collapse.  In  such  cases  no  amount  of  rezoning  can  prevent  mesh  tangling  and 
collapse  of  cells.  A  final  fault  of  the  rezone  solution  is  that  it  does  not  address  the  proper  questions.  A 
grid  deformation  must  force  a  less  accurate  approximation.  Rezone  algorithms  generally  seek  to 
preserve  a  reasonable  appearance  of  the  grid  and  may  obscure  the  question  of  what  is  the  best  approxi¬ 
mation  possible  given  the  current  Lagrangian  vertex  positions.  That  is,  under  the  guise  of  preventing 
inaccuracies  in  the  approximations,  they  may  become  the  very  vehicle  for  preserving  them. 

Although  these  remarks  have  been  cast  in  the  framework  of  finite  difference  techniques,  they 
apply  equally  to  Lagrangian  finite  element  techniques  (Strang  and  Fix,  1973).  Regardless  of  the  method 
used  to  calculate  vertex  velocities  and  accelerations,  a  Lagrangian  mesh  must  distort  and  if  a  rezone 
phase  is  applied  the  same  objections  as  to  its  accuracy,  adequacy  and  appropriateness  hold. 

Reconnecting 

The  solution  that  rezoning  offers  to  grid  deformation  is  vertex  motion  to  rectify  the  distortions. 
An  alternate  solution  is  illustrated  in  Figure  A3.  A  section  of  a  quadrilateral  mesh  about  a  shear  layer 
is  shown  in  Figure  A3a.  A  Lagrangian  calculation  quickly  leads  to  the  mesh  shown  in  Figure  A3b,  in 
which  mesh  connections  about  the  layer  no  longer  join  neighboring  vertices.  In  this  figure  one  grid  line 
has  been  reconnected  to  show  the  connection  which  is  now  appropriate.  For  a  periodic  system  all 
stretched  grid  lines  could  be  reconnected,  thereby  restoring  the  mesh  to  is  original  configuration  without 
moving  any  vertices.  In  general,  all  such  reconnections  either  are  inappropriate  or  cannot  be  made  due 
to  boundaries,  so  that  one  triangular  and  one  pentagonal  cell  remain.  Therefore  reconnection  on  a  qua¬ 
drilateral  grid  cannot  by  itself  resolve  the  problem  of  distorted  grids.  On  a  triangular  mesh,  however, 
there  are  no  such  complications.  As  shown  in  Figure  A4,  a  reconnected  grid  line  on  a  triangular  mesh 
still  results  in  two  triangular  cells.  This  technique  was  first  used  in  a  computer  code  by  Crowley  (1971) 
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and  represents  a  very  attractive  alternative  to  rezoning  for  triangular  grids.  There  is  no  Eulerian  vertex 
motion.  For  a  given  reconnection  only  one  grid  line  or  two  cells  are  affected,  instead  of  all  the  grid 
lines  and  cells  about  a  given  vertex  in  the  rezone  technique.  Therefore  at  most  diffusion  is  effectively 
over  just  one  cell  instead  of  two  in  either  direction.  The  simplicity  of  the  technique  is  particularly 
suited  to  considerations  of  how  best  to  construct  approximations  given  the  current  grid  positions.  As 
shown  below,  such  questions  result  in  very  attractive  conservative  algorithms  for  the  reconnection 
operation.  Its  major  similarity  with  rezoning  is  that  it  cannot  by  itself  solve  the  problem  of  fluids  evolv¬ 
ing  into  multiply-connected  regions  or  of  flows  at  stagnation  points.  However,  it  does  aid  in  the 
remedy.  The  number  of  grid  lines  meeting  at  any  vertex  can  be  reduced  to  three  by  reconnections, 
with  those  three  neighboring  vertices  forming  a  triangle  which  surrounds  only  that  vertex.  If  the  fluid 
is  accumulating  vertices  locally,  then  that  central  vertex  can  be  eliminated  with  the  three  grid  lines, 
leaving  only  the  surrounding  triangle.  The  result  is  the  desired  decrease  in  resolution  and  the 
avoidance  of  the  formation  of  thin,  narrow  triangles  near  the  point  of  converging  flow.  Conversely,  a 
vertex  may  be  added  inside  any  triangle  or  along  any  line  simply  by  providing  the  necessary  grid  lines 
to  other  vertices  within  the  affected  triangles.  Subsequent  reconnections  will  link  the  added  vertices  to 
their  neighbors.  In  this  way  the  transition  to  multiply  connected  regions  and  the  flow  near  stagnation 
points  can  be  handled  smoothly  merely  by  decreasing  or  increasing  resolution  where  appropriate.  The 
combination  of  grid  line  reconnection  with  vertex  addition  and  deletion  therefore  provide  a  means  of 
smoothly  restructuring  the  grid  without  recourse  to  Eulerian  vertex  movement. 

The  price  that  is  paid  for  this  flexibility  is  the  loss  of  any  global  ordering  for  the  vertices.  A 
reconnecting  grid  has  no  fixed  indexing.  Partly  because  of  this  difficulty  general  triangular  grids  have 
not  received  the  attention  given  to  regular  grids  and  are  not  as  well  understood.  Solution  procedures 
for  irregular  grids  are  currently  more  costly  and  the  range  of  available  numerical  techniques  is  much 
more  restricted.  Fortunately  grid  alterations  can  be  accommodated  in  finite-difference  schemes  with 
relatively  mild  calculational  repercussions.  For  finite  element  schemes  a  reassembly  of  the  system 
matrix  is  required  whenever  the  mesh  is  restructured.  For  that  reason  we  will  focus  on  finite  difference 
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techniques  from  this  point  on.  The  relative  lack  of  experience  in  differencing  over  general  triangular 
grids  has  historically  led  to  some  mistaken  impressions,  two  of  which  are  important  enough  to  address 
directly  —  their  "stiffness"  and  lack  of  accuracy. 

Placement  of  variables 

The  assignation  of  "stiff  to  triangular  grids  seems  to  have  evolved  from  an  improper  placement  of 
physical  variables  on  the  mesh.  In  a  quadrilateral  mesh,  there  is  a  one-to-one  correspondence  between 
cells  and  vertices.  For  example,  each  quadrilateral  cell  may  be  associated  with  its  upper  right  vertex. 
The  vertices  remaining  lie  on  boundaries,  and  their  variables  are  specified  by  boundary  conditions.  If 
that  same  quadrilateral  mesh  is  subdivided  into  a  triangular  mesh  by  drawing  diagonals  through  each 
cell,  it  is  clear  the  correspondence  between  cells  and  vertices  is  destroyed.  The  number  of  vertices  is 
the  same,  but  there  are  twice  the  number  of  cells.  On  the  quadrilateral  mesh  it  is  attractive  to  assign 
pressures  as  cell-centered  quantities,  since  pressure  forces  are  easily  calculated.  Such  a  positioning  on  a 
triangular  mesh  can  be  disastrous.  Twice  as  many  pressures  must  be  specified  than  in  the  well- 
formulated,  quadrilateral  case.  Numerically,  an  iterative  solution  for  pressures  would  converge 
extremely  slowly,  and  in  this  sense  the  system  of  differenced  equations  would  be  stiff.  If  the  idea  of 
conserving  cell  size  for  incompressible  flows  is  also  carried  over  to  triangular  grids,  the  algorithms  are 
stiff  in  a  much  worse  sense.  A  vortex  is  shown  in  Figure  A5,  where  the  arrows  designate  local  fluid 
velocities.  All  velocities  are  about  the  vortex  center,  and  there  exists  a  maximum  speed  at  some  radius 
from  the  center.  The  triangle  that  is  shown  will  therefore  eventually  invert,  since  the  vertex  with  the 
maximum  speed  must  pass  between  the  two  more  slowly  moving  vertices.  If  triangle  areas  are  con¬ 
served,  a  pressure  must  build  up  to  resist  this  tendency  until  the  speeds  of  all  vertices  are  commen¬ 
surate  with  cell  area  conservation.  The  presence  of  such  a  pressure  component  at  any  stage  of  the  cal¬ 
culation  is  totally  non-physical,  since  it  arises  from  improper  placement  of  variables  and  is  directly 
opposed  to  the  physical  flow.  Such  behavior  is  definitely  "stiff”  since  the  solution  will  tend  to  resemble 
solid  body  rotation  rather  than  vortex  motion.  Such  "stiffness"  of  triangular  grids  seems  to  be  due 
entirely  to  such  improper  differencing. 
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The  use  of  pressures  defined  on  vertices  will  illustrate  the  second  question,  that  of  accuracy  of 
difference  schemes  over  a  general  triangular  mesh.  First  let  us  incorporate  our  rather  novel  placement 
of  pressures  in  the  one-dimensional  case.  There  are  both  forward  and  backward  expansions  of  the  pres¬ 
sure  about  the  point  r, 

Pi+'  “  *  +  9%  Ax  +  T  0  A*2  +  A*3  +  0(Ax4)  (A1) 

and 

*-•  ■  p‘-  u,  4*  + 1 0  4*’ -  {0  4,3  +  0<Ax,) ■  <A2) 

Of  course  the  forward  and  backward  difference  approximations  can  be  obtained  directly  from  Equations 
A1  and  A2; 


and 


8p  ^  Pi+ 1  ~Pi 
f>Xj  Ax 


(A3) 


8p  =  Pi  ~  Pi-  i 
fix,  Ax 


(A4) 


Both  are  only  first-order  accurate  and  can  be  viewed  as  cell-centered  quantities.  If  instead  Equations 
A1  and  A2  are  subtracted,  the  result  is 

Pi+i  ~  Pi-\  “  2-&-  A x  +  y  J^T  Ax3  +  0(A*5).  (A5) 

The  centered  difference  approximation  is  then 


IE. 

8xt 


Pl+l  ~  Pi- 1 
2Ax 


(A6) 


which  is  second  order  accurate  and  vertex  centered.  For  the  case  of  variable  mesh  spacing,  we  can 


replace  Ax  by  Ax'  in  Equation  A2,  so  that  Equation  A5  becomes  instead 


A+i  -  pt- 1  -  -f^-  (Ax  +  Ax')  +  4  -—4-  (Ax2  -  Ax'2)  +  0( Ax3). 

OXf  2  Qa 

Therefore 

8p  _  Pi+ 1  ~  Pi-\ 

8xt  (Ax  +  Ax') 


(A7) 


(A8) 
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which  now  has  an  accuracy  between  first  and  second  order  depending  on  the  ratio  Ax/ Ax',  and  which  is 
evaluated  neither  directly  at  a  vertex  nor  the  cell  center.  If  we  rewrite  Equation  A4  using  primes  also, 
then  Equation  A8  can  be  rewritten  as 


-^-Ax,  +  Ax/ 
8 p  _  8x,  8x, 

8  X/  Ax,  +  Ax/ 


or 


Z 

8  p  _  <-i  oxi 

Sx,  IAx, 

That  is,  the  central  difference  can  be  obtained  by  an  "area"  weighted  sum  of  the  forward  and  backward 


(A9) 


differences. 


This  result  carries  over  directly  to  general  triangular  grids  in  two  dimensions.  We  can  define  the  pres¬ 
sure  gradient  there  as 

(A10) 

/-I  2A 

where  Vp  is  the  first  order  accurate,  finite  difference  approximation  to  the  gradient  which  is  evaluated 
at  the  triangle  centroid,  A  is  the  triangle  area,  z  is  a  unit  vector  in  the  direction  of  the  neglected  coordi¬ 
nate  and  the  sum  extends  over  the  three  triangle  vertices.  The  analogue  of  Equation  A9  is 

t^PjAj 

Vpf  ~  (AIl) 

which  for  special  geometries  is  second  order  accurate  or  higher.  In  general  it  is  less  then  second  order 
accurate,  but  performs  reasonably  close  to  second  order  accuracy  for  a  general  mesh  provided  that  the 
triangle  areas  are  roughly  equal.  In  that  case  the  error  is  determined  by  a  formula  similar  to  Equation 
A7,  an  algebraic  sum  of  squared  triangle  altitudes.  The  worst  that  can  be  achieved  is  first  order  accu¬ 
racy,  which  is  obtained  only  in  the  degenerate  case  of  a  zero  area  triangle  at  the  vertex.  This  implies 
that  care  must  be  taken  in  evaluating  boundary  conditions,  but  this  task  is  alleviated  by  having  the  pres¬ 
sures  specified  at  the  boundaries.  In  a  cell  centered  scheme,  the  pressures  are  located  half  a  cell  away. 


i 
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and  boundary  conditions,  particularly  at  free  surfaces,  are  more  difficult  to  implement.  Accuracy  in  the 
interior  of  the  fluid  is  therefore  diminished  primarily  by  narrow  triangles.  As  shown  below,  this  restric¬ 
tion  is  not  too  serious  for  a  reconnecting  grid  since  the  grid  can  be  made  to  reconnect  to  preserve  regu¬ 
lar  triangles.  In  cases  where  this  is  not  possible  (near  interfaces,  for  example),  the  grid  addition  or 
deletion  vertices  can  be  used  to  regularize  the  mesh. 

The  conclusion  is  that  the  "stiffness/'  and  ”  diminished"  accuracy  of  approximation  on  the  triangular 
grid  are  both  in  large  measure  due  to  improper  placement  of  physical  variables  on  the  mesh.  Pressures 
as  cell-centered  quantities  force  not  only  a  stiff  solution,  but  make  second-order  accurate  approxima¬ 
tions  difficult  to  achieve.  Two  other  comments  should  be  made  about  Equation  All.  The  first  is  that  it 
can  be  used  to  recover  all  the  usual  gradient  approximations  for  regular  grids.  An  extension  of  the 
definition  by  taking  the  divergence  of  both  sides  of  Equation  All  also  yields  all  the  usual  Laplacian 
approximations  (see  below).  The  second  feature  of  the  equation  is  that  the  vertex  pressure  gradients 
can  be  viewed  either  as  a  sum  over  triangular  gradients  or  in  the  more  conventional  way  of  a  vertex 
sum.  The  former  offers  the  possibility  of  vectorizing  the  equations  even  though  the  mesh  has  no 
inherent  order,  a  fact  which  greatly  enhances  the  utility  of  the  method  for  the  new  vector  computers. 

Numerical  Algorithms 

The  computer  code  SPLISH  is  a  two-dimensional  Lagrangian  fluid  dynamics  code  for  incompressi¬ 
ble  fluids  which  was  developed  in  accordance  with  the  philosophy  of  the  previous  section  (Boris  et.  al., 
1975,  and  Fritts  and  Boris,  1977).  The  basic  equations  for  incompressible  fluid  flow  are 

p-^--Vp-pgP  (A12) 

at 

and 

V  •  V  -  0  (A13) 

where  the  fluid  density  p,  pressure  p  and  velocity  V  are  assumed  to  vary  only  with  x  and  y.  With  pres¬ 
sures  specified  at  the  vertices,  Vp  is  evaluated  over  triangles,  and  Equation  A12  can  easily  be  updated 
implicitly  or  explicitly  if  velocities  are  considered  to  be  triangle-centered.  This  placement  of  velocities 
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as  cell  quantities  and  pressures  at  vertices  is  apparently  unique  to  SPLISH  and  is  the  direct  opposite  of 
the  usual  placement.  In  what  follows  the  subscript  /'  will  denote  a  vertex-centered  quantity  and  j  a 
triangle-centered  quantity.  In  SPLISH  the  integration  of  velocities  uses  a  split  step  algorithm  whereby 
the  velocities  are  advanced  one  half  timestep  (Equation  A14),  the  grid  is  advanced  a  full  timestep 
(Equation  A 16),  and  then  the  velocities  advanced  forward  the  other  half  timestep  (Equation  A18). 


(A14) 

v?-j<vr+viO. 

(A15) 

1 

X"  X,0  +  St  V,2 , 

(A16) 

v/  =  R({Xf},  (X/1))  ■  \j, 

(A17) 

(A  18) 

The  vertex  velocity  V"  appearing  in  Equation  A15  is  obtained  from  the  area-weighted  V"  from  the  pre¬ 
vious  iteration, 

IVM 

V,”-  V.  (A19) 

The  advantage  of  using  triangle  centered  velocities  is  the  ease  in  conceptualizing  and  expressing  conser¬ 
vation  laws.  Because  of  the  paucity  of  experience  in  formulating  algorithms  over  a  general  triangular 
grid,  we  employed  the  control  volume  approach,  which  uses  an  integral  formulation  to  derive  the 
difference  algorithms.  Equation  A17  is  the  first  manifestation  of  this  approach.  It  reflects  numerically 
the  fact  that  the  triangle  velocities  must  rotate  and  stretch  as  the  grid  rotates  and  stretches.  The 
transformation  R  is  derived  by  considering  the  circulation  about  each  vertex.  The  boundaries  of  a  ver¬ 
tex  cell  can  be  defined  by  apportioning  each  triangle  equally  to  each  of  its  vertices.  One  way  of  achiev¬ 
ing  this  partitioning  is  to  draw  bisectors  from  the  mid  point  of  each  side  to  the  opposite  vertex,  as 
shown  in  Figure  A6.  The  side  bisectors  meet  at  the  centroid  of  the  triangle  and  divide  its  area  into  six 
smaller  equal  area  triangles.  The  areas  of  the  two  small  triangles  adjacent  to  each  vertex  are  assigned  to 
that  vertex.  The  vertex  cell  of  Figure  A6  is  constructed  by  summing  over  all  the  surrounding  triangles. 
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Therefore  the  area  of  a  vertex  cell  may  be  defined  as 

Ac  -  £  \  Aj  (A20) 

j  3 

where  the  sum  extends  over  all  adjacent  triangles,  exactly  as  in  Equation  All.  With  this  definition 
Equation  All  and  Equation  A19  become 

VpjAj 

Vft - J— -  (A21) 

Ac, 

and 

VjAj 

y.  =  ±J- - .  (A22) 

AC/ 

Since  the  triangle  velocities  are  constant  over  the  triangle,  the  circulation  taken  about  the  boundary  of 

the  vertex  cell  is  straightforward.  Circulation  is  conserved  about  each  of  the  three  triangle  vertices  by 

the  transformation  R  of  Equation  A17.  For  constant  velocities  the  transformation  is  linear  and  is  given 
by  the  three  equations; 

v/  •  (Xf-  Xf)  =  v/  •  (Xf-  Xf), 

V/  •  (Xf  -  Xf)  =  V/  •  (Xf  -  Xf),  (A23) 

v/  •  (Xf-Xf)  =  V/  •  (Xf-Xf). 

where  the  subscripts  1,  2  and  3  refer  to  the  vertices  associated  with  triangle  j.  This  transformation 
ensures  that  the  vorticity  integral  calculated  about  any  interior  vertex  is  invariant  under  the  advance¬ 
ment  of  the  grid.  It  is  easy  to  show  that  the  Vp  and  gravity  terms  cannot  alter  the  vorticity  either  since 

V  pi 

numerically  VxVp  =  0  and  gravity  is  a  constant.  Only  the  term  can  change  vorticity,  exactly  as 

dictated  by  the  physics.  Since  the  transformation  R  is  time  reversible,  so  are  Equations  A14-A18,  so 
that  the  entire  algorithm  advances  vertex  positions  and  velocities  reversibly  while  evolving  the  correct 
vorticity  about  every  interior  vertex.  No  numerical  generation  of  vorticity  can  occur,  a  rather  unique 
feature  of  SPLISH  among  Lagrangian  codes.  The  advantage  of  calculating  triangle  velocities  is  clearly 
evident  in  the  simplicity  of  the  transformation  R,  as  well  as  in  the  reconnection  algorithms,  as  shown 
below. 
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The  pressure  p"  in  Equation  A18  are  derived  from  the  condition  that  the  new  velocities  V]  should 
be  divergence  free  at  the  new  timestep,  satisfying  Equation  A13.  The  pressure  Poisson  equation  is 
derived  from  Equation  A18  by  setting  (V  •  Vf),  =  0,  to  obtain  pressure  p"  such  that 

i£(V  •  -L  (Vp)/),  -  (V  •  (A24) 

L  Pj 

The  right  hand  side  of  Equation  (A24)  is  the  numerical  analogue  of  the  V  •  (V  •  VV)8 1  term  which 
arises  when  the  divergence  of  Equation  A12  is  taken.  Both  terms  in  Equation  A24  are  simple  to  evalu¬ 
ate  since  the  divergence  is  taken  over  triangle  centered  quantities.  The  paths  are  the  "surfaced  bound¬ 
ing  the  vertex  volume  of  Figure  A6,  where  the  normal  is  directed  outward  from  the  vertex.  The  Pois¬ 
son  equation  (Equation  A24)  that  resulted  from  this  integration  has  two  advantages.  First  it  is  derived 
from  V2<f>  =  V  •  V<£  as  in  the  continuum  case.  Secondly  the  left-hand  side  results  in  the  more  fami¬ 
liar  second  order  accurate  templates  (such  as  the  five-point  formula)  for  the  Laplacians  for  homogene¬ 
ous  fluids  and  regular  mesh  geometries. 

In  summary  the  finite  difference  formulas  for  SPLISH  are  derived  using  a  control  volume 
approach.  Specifications  of  pressure  at  vertices  leads  naturally  to  the  choice  of  positioning  velocities  at 
triangles.  Although  pressure  gradients  are  constant  over  triangles,  the  resultant  algorithms  are  expected 
to  be  nearly  second  order  accurate  since  vertex  velocities  are  derived  from  pressure  gradient  forces 
through  sums  about  vertices,  which  in  effect  takes  the  central  differences.  The  code  has  been  tested 
extensively  on  finite  amplitude  standing  waves  and  has  been  shown  to  be  basically  second  order  accu¬ 
rate  by  the  variation  in  period  with  mesh  size  (Fritts,  1976a  and  Fritts  and  Boris,  1979). 

Grid  restructing  algorithms 

The  derivation  of  the  reconnection  and  vertex  addition  and  deletion  algorithms  is  accomplished  in 
the  same  manner  through  the  control  volume  approach  and  the  use  of  triangle  velocities.  The  accuracy 
of  a  computer  code  which  uses  a  reconnecting  grid  is  determined  by  two  aspects  of  the  algorithms— how 
accurately  post-reconnection  physical  variables  are  chosen  and  when  the  reconnections  occur.  As  men¬ 
tioned  above,  a  reconnection  offers  the  possibility  of  much  less  diffusion  since  fluid  passes  through  only 
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one  grid  line  for  reconnections,  whereas  four  grid  lines  are  involved  in  every  vertex  rezone  movement. 
Every  mesh  line  can  be  viewed  as  one  diagonal  of  the  quadrilateral  formed  by  the  triangles  to  either 
side.  A  reconnection  merely  chooses  the  other  diagonal.  During  a  reconnection  the  smallest  definable 
cell  is  the  quadrilateral,  not  the  two  triangles,  and  it  is  necessary  to  ensure  that  quadrilateral  properties 
are  unchanged  during  a  reconnection.  That  is,  the  quadrilateral  is  a  control  volume  over  which  certain 
physical  variables  are  conserved.  As  shown  in  Figure  A7,  a  reconnection  alters  the  vertex  cells  for  each 
of  the  four  quadrilateral  vertices.  To  keep  vorticity  and  divergence  conserved  the  portions  of  the 
integrals  J*  V  •  \dV  and  J*  V  x  V  •  d\  about  each  vertex  which  lie  within  the  quadrilateral  must  be 
the  same  before  and  after  the  reconnection.  Since  there  are  four  velocity  components  to  be  specified 
after  the  reconnection,  it  would  seem  possible  only  to  keep  either  the  divergence  or  the  vorticity  about 
each  vertex  constant,  but  not  both.  However,  the  divergence  and  vorticity  equations  are  not  all 
independent  and  it  is  possible  to  exactly  conserve  both  vorticity  and  divergence  about  each  vertex  with 
a  unique  set  of  triangle  velocities.  This  same  solution  also  conserves  the  quadrilateral  velocity  and  has 
the  added  feature  of  time  reversibility.  Re-reconnecting  a  line  yields  the  identical  grid  and  physical 
variables  both  on  vertices  and  triangles.  This  is  extremely  desirable  since  the  basic  finite-difference 
equations  were  also  time  reversible,  as  are  the  physical  equations  themselves. 

The  question  of  when  to  reconnect  remains  unresolved.  As  mentioned  above,  the  accuracy  of  a 
general  triangular  mesh  is  diminished  by  narrow  triangles.  With  reconnections,  accuracy  can  be 
recovered  by  ensuring  that  small  angles  are  preferentially  eliminated.  There  are  many  ways  of  quantify¬ 
ing  such  an  algorithm.  The  one  we  have  chosen  arises  naturally  from  the  pressure  Poisson  equation. 
Since  the  equation  is  solved  in  SPLISH  by  iteration,  it  is  desirable  that  the  convergence  of  the  iteration 
be  as  rapid  as  possible.  Mathematically,  convergence  is  assured  if  the  equation  exhibits  diagonal  domi¬ 
nance.  For  a  general  triangle  mesh  diagonal  dominance  is  obtained  only  if  all  coefficients 

a  -  1/2  (cot  9~  +  cot  9+ )  (A25) 

are  positive,  where  9+  and  9~  are  defined  in  Figure  A8a.  For  positive  area  triangles  9+  and  9~  are  both 

between  0°  and  180°,  so  that  each  term  is  negative  only  when  9*  +  9~  >  180°,  since 

IS 


MINER,  GRIFFIN,  RAMBERG,  AND  FRITTS 


sin(0+  +  9  ) 


(A26) 


2sin0+sin0— 

In  SPLISH  the  reconnection  algorithm  is  chosen  to  be  exactly  the  algorithm  needed  to  preserve  diago¬ 
nal  dominance.  If  9+  +  9~  is  greater  than  180°,  the  grid  line  is  reconnected  as  shown  in  Figure  A8b. 
The  new  angles  0  +  and  9  ~  must  seem  to  less  than  180°  since  9+  +  9~  +  9  +  +  9  ~  is  the  sum  of  the 
interior  quadrilateral  angles,  which  must  be  360°.  Therefore  the  reconnection  algorithm  is  unique  as 
well  as  straightforward.  It  also  preferentially  eliminates  small  angles  for  triangles,  si®ce  the  diagonal  is 

chosen  which  divides  the  largest  opposing  angles. 


The  second-order  accuracy  of  SPLISH  may  be  expected  to  be  preserved,  because  the  reconnection 
algorithms  are  specified  to  ensure  diagonal  dominance  and  eliminate  small  angle  triangles.  As  yet  no 
test  has  been  made  of  the  accuracy  of  the  reconnection  algorithms  for  the  complete  range  of  gridding 
situations.  This  is  mainly  because  reconnections  cannot  themselves  ensure  second-order  accuracy, 
since  flows  near  stagnation  points  must  force  narrow  triangles  regardless.  That  is,  vertex  addition  and 
deletion  must  be  employed  at  the  same  time.  The  reconnection  algorithms  have  been  tested  exten¬ 
sively,  however,  through  the  Kelvin-Helmholtz  instability.  In  the  linear,  non-linear  and  turbulent 
regimes,  the  algorithms  have  been  shown  to  provide  accurate  calculations  by  comparisons  with  both 
theory  and  experiment  (Fritts,  1976b  and  Fritts  and  Boris,  1979). 

The  derivation  of  the  remaining  grid  restructuring  algorithms  proceeds  in  exactly  the  same 
manner  as  above.  For  vertex  addition  satisfaction  of  conservation  integrals  is  particularly  simple.  A 
vertex  added  at  the  centroid  of  a  triangle  subdivides  that  triangle  into  three  smaller  triangles.  If  the 
new  triangle  velocities  are  all  the  same  as  the  velocity  of  the  subdivided  triangle,  all  conservation  laws 
are  trivially  satisfied.  Since  the  reconnection  algorithm  is  also  conservative,  subsequent  reconnections 
to  other  vertices  ensure  that  the  only  effect  of  the  addition  is  an  increase  in  resolution. 


The  case  is  not  as  obvious  for  vertex  deletion.  Reconnections  can  be  used  to  surround  any  inte¬ 
rior  vertex  within  a  triangle.  It  can  be  shown,  however,  that  once  the  vertex  is  surrounded  by  a  trian¬ 
gle,  the  motion  of  that  triangle  is  not  altered  if  the  vertex  is  removed  and  the  new  larger  triangle  given 
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a  velocity  which  is  the  area-weighted  sum  of  the  old  velocities. 


AJ/t  —  A  iVt  +  A  2V2  +  A  3V3.  (A27) 

Such  a  substitution  also  conserves  vorticity  exactly,  and  effects  a  redistribution  in  accordance  with  area 
coordinates.  Figure  A9  illustrates  the  triangles  before  and  after  vertex  removal.  If  £4  is  the  vorticity 
about  vertex  4  before  removal,  then  the  vorticity  about  each  of  the  other  three  vertices  is  increased  by 
an  amount  f where 

f'l  “  AjtjA, 

f  j  -  A&JA,  (A28) 

fa  =  A&JA, 

where  +  £'2  +  £'3  *=  £4  since  A,  +  Aj  +  Ak  =  Ah  Therefore  total  vorticity  is  conserved  and  redistri¬ 
buted  in  a  reasonable  and  natural  manner.  Since  the  behavior  of  the  divergence  is  governed  by  a  simi¬ 
lar  set  of  equations  and  conservation  of  momentum  is  ensured  by  Equations  A27,  the  conservation  of 
the  flow  variables  is  guaranteed  and  a  loss  in  resolution  is  the  primary  effect  of  deletion. 

The  vertex  addition  and  deletion  routines,  together  with  the  reconnection  algorithms  are  incorporated 
in  SPLISH  through  driving  routines  which  test  for  large  or  small  grid  lines  and  large,  small  or  skewed 
triangles,  as  well  as  special  tests  at  boundaries  and  interfaces.  The  resultant  code  automatically  restruc¬ 
tures  the  grid  under  constraints  in  the  form  of  a  maximum  allowable  skewness  and  maximum  and 
minimum  areas  and  line  lengths.  The  complete  routine  has  been  tested  using  the  Rayleigh-Taylor  ins¬ 
tability  and  in  flows  about  hydrofoils.  Preliminary  results  of  these  tests  are  extremely  satisfying,  pro¬ 
viding  the  first  totally  Lagrangian  calculations  of  these  phenomena  at  long  times  (Fritts,  1977  and 
1979). 

The  tests  performed  using  this  code  for  the  wave-structure  interaction  program  are  specifically  aimed  at 
the  nonlinear  free  surface  effects  to  answer  the  question  of  whether  the  code  does  indeed  still  behave 
as  a  basically  second-order  code  and  if  the  boundary  conditions  used  thus  far  are  realistic. 
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Figure  Al.  (a)  Grid  connections  for  a  regular  mesh,  and  (b)  a  distorted  mesh 
formed  by  shear  motion  of  the  central  column  relative  to  its  neighbors. 


(b) 


•  •  • 


Figure  A2.  (a)  Vertices  used  in  higher  order  approximations,  (b)  The  vertex 
positions  on  the  mesh  distorted  by  shear. 
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(a)  (b) 


Figure  A3,  (a)  A  quadrilateral  mesh  about  a  shear  inter¬ 
face.  (b)  Distortions  induced  in  the  Lagrangian  mesh  by 
the  shear.  A  reconnection  has  been  made  to  join  nearest 
neighbors  for  a  single  vertex. 


Figure  A4.  (a)  A  triangular  grid  formed  by  drawing  diagonals 
through  the  mesh  of  Figure  A3a.  (b)  The  reconnection  of  Fig¬ 
ure  A3a  on  a  triangular  mesh. 
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Figure  A5.  An  isolated  vortex  about  a  central  vertex.  Velocities  of  neighboring  vertices  reach  a  maximum  at 
some  radius  from  the  central  vertex  and  diminish  further  out. 
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V3 


Figure  A6.  A  definition  of  a  vertex  cell  for  a 
general  triangular  mesh. 
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Figure  A8.  (a)  Definition  of  the  angles  9*  and  9~  for  the  diagonal  line 
drawn  from  j  to  /.  (b)  The  angles  @'+  and  O'-  formed  by  connecting  the 
other  quadrilateral  diagonal. 
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Appendix  B 

WAVE  INDUCED  FLUID  PARTICLE  MOTIONS 

It  is  straightforward  to  compare  the  Lagrangian  numerical  computations  with  classical  solutions  to 
the  equations  of  motion  for  inviscid,  irrotational  flow.  Moreover,  as  demonstrated  originally  by  Stokes 
(1847)  and  more  recently  by  Phillips  (1966),  Newman  (1977)  and  Longuet-Higgins  (1980),  it  is  con¬ 
venient  to  derive  the  properties  of  the  wave  motion  from  existing  Eulerian  solutions  and  their  Lagran¬ 
gian  analogues. 

A  two-dimensional  coordinate  system  is  chosen  with  y  —  0  in  the  plane  of  the  undisturbed  free 
surface  and  with  the  y  axis  positive  upwards.  The  positive  x  axis  is  taken  from  left  to  right  in  the  direc¬ 
tion  of  wave  motion.  If  the  velocity  of  a  fluid  element  in  Lagrangian  terms  is  denoted- by  the  vector 
u/(a,  /),  with  the  initial  condition  x  =  a  at  /  —  0,  then  the  position  of  the  element  at  a  later  time  is 

x  -  a  +  J0  U/ (a,  t)  dt.  (Bl) 

The  velocity  u/  is  given  by  the  Taylor  series  expansion 

U/(a,  t)  -  u(a,  t)  +  |j*0  u/(a,  t)  rflj  ■  Vau(a,  t)  +  ,  (B2) 

where  «(a,  /)  is  the  Eulerian  vector  velocity  and  V„  is  the  gradient  operator  in  the  direction  a.  To  first 
order  Stokes’  solutions  for  the  velocity  fields  are  identical  in  Lagrangian  and  Eulerian  terms,  so  that 

u/(a,  t)  -  u(a,  t )  +  |j*0  n(a,  t )  d (J  •  Va  u(a,  t),  (B3) 

correct  to  second  order  in  the  surface  displacement  amplitude.  A  wave  with  the  displacement  profile 
(where  H  is  the  wake  amplitude) 

Tj  -  H  cosikx  -  (Tt)  (B4) 

from  the  undisturbed  free  surface  is  satisfied  by  the  first-order  Eulerian  velocity  potential 

»(*■>)-  <7"C”h  +  J>  sintfac  -  <r,).  <SS) 

k  sinh  kd 

and  the  corresponding  dispersion  relation 
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<r2  -  gk  tanh(Ad)  (B6) 

for  water  of  uniform,  finite  depth  d.  Here  a  denotes  the  wave  frequency  <r  -  In!  T,  where  T  is  the 

wave  period;  and  k  denotes  the  wave  number  k  -  2tt/A,  where  A.  is  the  wavelength.  For  deep  water 

waves  with  d/k  >  0.5  the  dispersion  relation  reduces  to  a-2  -  gk.  The  expressions  for  the  velocity 

potential  and  dispersion  relation  satisfy  the  continuity  equation  V  •  u  •  0,  which  reduces  to  Laplace’s 

equation 

V2</>  -  0  (B7) 

since  u  —  and  v  =  Also  satisfied  are  the  linearized  boundary  conditions 
dx  9y 

Kinematic:  (B8) 

9y  at 

Dynamic:  i)  —  — -  (B9) 

g  9 1 

at  the  undisturbed  free  surface,  y  —  0.  The  normal  velocity  is  assumed  to  vanish  at  the  plane  bottom, 
thus  satisfying  the  additional  boundary  condition 

v  -  M  -  0  at  y  =  —d.  (BIO) 

9y 

The  Eulerian  fluid  velocities  which  are  obtained  from  this  system  of  equations  are 

d<f>  <tH  cosh  k(y  +  d)  ,, 

u  =  - - —■  - cos  (toe  -  a-t)  (Bll) 

ox  stnh  kd 

and 


d4>  o-H  sinh  k(y 
9y  “ 


In  Cartesian  form  the  vector  equation  u(  — 
the  initial  position  a  -  0c„,  y0)  is  given  by 


cosh  kd 


dx  dy 
dt’  dt 


—  —  sin(Ax  -  err).  (B12) 

for  the  Lagrangian  velocity  of  a  fluid  particle  with 


-rr  -  u{x0,  y0,  t)  +  (x  -  Xo )  +  (y  -  y„)  (B13) 

”  v^0-  y0,  t)  +  (x  -  Xo)  +  (y  -  y0)  -J-^,  (B14) 

correct  to  second  order.  Here  the  velocities  u  and  v  are  given  by  Equations  Bll  and  B12  evaluated  at 
the  point  (x0,  y„),  and  the  first-order  trajectories  are  given  by 
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and 


X  -  x0 


y0.  t )  dt 


(B15) 


y  -  y0  “  J0  v(x0,  y0,  t)  dt.  (B16) 

An  interesting  finding  that  is  derived  from  Equations  B13  and  B14  is  the  existence  of  a  second-order 
mean  motion  of  the  fluid  particles  in  the  x  direction.  The  particle  paths  are  not  closed  and  elliptical  as 
in  the  first-order  solution,  but  instead  there  is  a  slow,  depth-dependent  drift  of  the  fluid  particles  in  the 
direction  of  the  wave  motion.  There  is  no  mean  flow  in  the  y  direction  and  v  remains  periodic,  so  that 
one  obtains  the  mean  velocities 


crkH1  cosh  2k  (y0  +  d) 
u,=  — 


2  sinh2  kd 

V/  —  0. 

This  second-order  mean  motion  in  the  direction  of  the  wave  propagation  is  often  called  the  "Stokes 


(B17) 

(B18) 


drift." 
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Appendix  C 

EXPERIMENTAL  EQUIPMENT  AND  METHODS 


The  wave  channel  employed  in  the  present  experiments  is  housed  in  Building  A59  at  NRL.  An 
overview  of  the  channel,  the  wind  generator,  and  the  mechanical  wave  generator  is  shown  in  Figure  Cl. 
A  photograph  of  the  bottom-sealed  half  cylinder  in  the  channel  is  shown  in  Figure  C2.  The  channel  is 
30.5m  (100  ft)  in  length  and  1.22m  (4  ft)  wide,  and  the  maximum  water  depth  is  1.07m  (3.5  ft).  The 
regular  waves  required  for  the  experimental  program  discussed  in  this  report  were  generated  by  a 
bulkhead-type  mechanical  wavemaker  driven  by  a  variable-speed  DC  motor  and  a  scotch  yoke  mechan¬ 
ism.  A  sloping  beach  is  located  in  the  channel  at  the  opposite  end  from  the  wavemaker  in  order  to 
absorb  the  incident  waves.  The  beach  is  an  open  frame  construction  covered  with  porous  acrylic  plates 
and  a  rubberized  hair  blanket.  This  combination  results  in  beach  reflection  coefficients  ( HR/Hn  of  0.05 
or  less  over  most  of  the  wavelength  range  of  interest  in  the  present  experiments.  Larger  beach 
reflections  were  encountered  in  the  longer  wavelength  experiments  and  this  point  is  discussed  later  in 
this  appendix.  Details  of  the  channel  and  the  laboratory  facility  are  given  in  another  NRL  report  (Cl). 


Significant  beach  reflections  (>  2%) ,  in  the  absence  of  a  submerged  object,  are  readily  detected  as 
a  modulation  in  waveheight  along  the  channel  (C2).  In  its  simplest  form,  this  amplitude  modulation 
arises  from  the  superposition  of  three  linear  progressive  wave  trains:  the  incident  waves,  the  partially 
reflected  waves,  and  the  essentially  total  reflection  of  the  latter  waves  by  the  wavemaker.  The  beach 
reflection  coefficient  Rb  in  this  case  is  given  by  the  modulation  extremes  arranged  as  (C2) 


R„- 


Hrmx  Hmin 
^max  4"  ^min 


(Cl) 


and  clearly  requires  that  any  additional  or  secondary  beach  reflections  are  negligible,  i.e.  Rb  «  1.  The 
same  procedure  can  be  used  to  determine  the  cylinder  reflection  coefficient  Rc,  however,  secondary 
reflections  tend  to  be  significant  depending,  of  course,  on  the  magnitude  of  R?.  In  fact  one  can  easily 
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envision  between  the  cylinder  and  the  wavemaker  a  progression  of  higher  order  standing  waves  whose 
relative  importance  grows  rapidly  with  increasing  Rc. 

A  further  complication  arose  in  the  present  experiments  because  the  beach  performance  was 
degraded  at  the  longer  wavelengths  and  lower  water  levels.  Under  those  conditions,  the  waves 
transmitted  by  the  cylinder  were  partially  reflected  by  the  beach  and  returned  to  strike  the  "back"  of  the 
submerged  half-cylinder.  The  initial  transmitted  waves,  denoted  by  the  transmission  coefficient  T  and 
satisfying  the  energy  balance  (C3) 

P  =  1  —  Rc2,  (C2) 

are  themselves  transmitted  and  reflected  by  the  cylinder.  This  gives  rise  to  additional  standing  wave 

components  in  front  of  and  behind  the  cylinder.  The  essential  difficulty  in  this  type  of  laboratory 
experiment  is  clear.  Each  wave  train  that  is  incident  upon  the  cylinder  produces  two  outgoing  wave 
trains  which  are  returned  either  partially  or  totally  to  establish  two  additional  standing  wave  components 
and  two  additional  incident  waves  for  the  next  interactions.  It  is  important  to  note  that  only  part  of  this 
difficulty  is  associated  with  large  beach  reflections  and,  in  fact,  large  corrections  will  be  required  for 
seemingly  simply  cases.  For  example,  consider  a  case  where  Rb  —  0.05,  acceptable  by  common  stan¬ 
dards,  and  Rc  =  0.10.  A  modulation  measurement  on  the  wavemaker  or  upstream  side  of  the  half 
cylinder  includes  not  only  a  standing  wave  due  to  cylinder  reflection  but  also  a  comparable  standing 
wave  due  to  transmission  past  the  cylinder  and  to  beach  reflection.  An  expression  for  the  observed 
modulation  coefficient  Rm  in  this  instance  can  be  found  by  an  extension  of  the  earlier  simple  method 
(C2)  and  is  of  the  form 

R2  -  Rc2  +  (R„  r2)2  +2 RcRb fW  (C3) 

where  <f>  is  the  phase  angle  between  the  two  standing  waves.  The  difference  between  the  observed 

modulation  coefficient  and  the  actual  cylinder  reflection  for  the  above  example  can  be  nearly  as  much 

as  ±  50  percent.  This  equation  was  employed  to  determine  the  cylinder  reflection  coefficients  shown  in 

Figure  32.  The  waveheight  modulations  and  relative  phases,  with  and  without  the  cylinder  in  place, 

were  obtained  from  stripchart  recordings  of  the  output  from  a  pair  of  traversing,  capacitance-type  wave 
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gages.  For  the  present  ranges  of  wave  steepnesses  and  reflection  coefficients,  second  harmonic  correc¬ 
tions  to  the  modulation  measurement  (C3)  were  considerably  smaller  than  the  above  correction  and 
therefore  were  neglected.  It  is  important  to  note,  however,  that  the  above  method  neglects  secondary 
reflections  by  the  assumptions  that  Rf,  Rc2,  and  Rt,Rc  are  all  much  less  than  unity.  This  was  not  the 
case  during  the  experiments,  and,  in  fact,  the  discrepancies  in  Figures  33  and  36  are  comparable  to 
these  quantities  for  each  case. 

The  wave-induced  pressures  on  the  submerged  half-cylinder  were  obtained  from  nineteen 
equally-spaced  =10°)  pressure  taps  located  around  the  circumference  of  the  cylinder  at  its  midsec¬ 
tion.  The  individual  taps  were  connected  to  a  differential  pressure  transducer  (Celesco  P7D)  by  a  24 
poskwm,  remotely-operated,  rotary  pressure  switch  (Scanivalve).  The  extra  switch  positions  were  used 
to  obtain  the  wave-induced  pressures  away  from  the  cylinder  and  to  purge  the  lines  of  air  by  means  of  a 
solenoid  valve  arrangement.  The  resolution  of  the  pressure  measurements  was  enhanced  by  approxi¬ 
mately  balancing  the  relatively  large  hydrostatic  pressure  with  a  separate  column  of  water  connected  to 
the  reference  side  of  the  pressure  transducer.  Calibration  of  the  pressure  sensor  was  achieved  by  vary¬ 
ing  this  reference  pressure  over  a  known  range.  The  pressure  signals  were  digitized  and  processed  for 
each  tap  position  by  means  of  a  Hewlett-Packard  5420A  Digital  Signal  Analyzer.  The  signal  magnitudes 
and  their  relative  phases  were  obtained  from  the  autospectra  and  cross-correlation  (phase)  functions, 
respectively.  The  overall  accuracy  of  the  measurement  system  was  better  than  5  percent.  Instantane¬ 
ous  pressure  distributions  were  reconstructed  from  the  measured  distributions  of  the  pressure  fluctua¬ 
tion  magnitudes  and  relative  phases.  The  presence  of  additional  waves  in  the  channel,  as  described  ear¬ 
lier,  often  had  a  large  influence  on  the  measured  pressures  and  this  was  taken  into  account  by  removing 
the  unwanted  wave  contribution(s)  according  to  linear  wave  theory  and  by  using  the  observed  wave 
modulation  envelopes  both  upstream  and  downstream  of  the  cylinder. 
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Figure  Cl  -  An  overview  of  the  NRL  wave  channel  and  laboratory. 
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Figure  C2  -  A  photograph  of  the  bottom  seated  half-cylinder  in  the  wave  channel. 
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